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CHAPTER  1 

INTRODUCTION  AND  SUMMARY 

1.1  INTRODUCTION 

This  document  constitutes  the  Final  Technical  Report  for  contract  F49620-90-C -00 1 6 
Application  of  Gabor  Representations  to  Military  Problems.  All  work  performed  by  Atlantic 
Aerospace  Electronics  Corporation  under  the  Final  year  of  this  contract  is  reported  herein. 

The  outline  of  this  report  is  as  follows.  The  remainder  of  this  chapter  contains  short 
summaries  of  the  work  items.  Chapter  2  reviews  theoretical  developments.  Chapter  3  discusses 
the  software  tasks  and  Chapter  details  the  applications  areas  investigated.  There  are  three  of 
these:  (4.1)  application  of  morphology  to  specific  emitter  analysis;  (4.2)  use  of  Gabor  transforms 
in  automatic  target  recognition  and  (4.3)  minimum  dimension  Gabor  representations. 

1.2  SUMMARY 

122.1  Theoretical  Developments 

Discussion  of  theory  is  presented  in  two  parts:  analytical  results;  and  algorithms. 

122.1.1  Theory 

Theory  topics  investigated  include  a  Gabor  sampling  theorem,  further  studies  on  accuracy  and 
stability  of  the  representation  and  relations  between  matrix  method  stability  and  signal  duration. 

122.1.1.1  Gabor  Sampling  Theorem 

In  section  2.1.1  we  show  a  sampling  theorem  applicable  to  Gabor  expansions  of  bandlimited 
functions,  derived  rather  directly  by  Shannon's  technique.  The  theorem  illuminates  relationships 
that  exist  between  the  number  of  degrees  of  freedom  in  the  signal  and  window  functions  and  the 
minimum  number  of  Gabor  coefficients  that  must  be  retained  to  have  a  faithful  expansion. 

1.2.1.12!  Accuracy  and  Stability 

In  section  2.1.2  we  report  on  the  stabilization  of  the  Gabor  transform  under  situations  where  a 
window  choice  desirable  from  an  application  point  of  view  leads  to  a  near  singular  mapping. 
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Methods  applicable  to  both  versions  of  the  transform  used  in  GS PS  CZak  and  matrix)  are 
presented.  In  addition,  some  thought  has  been  given  to  the  use  of  discrete  representations  in 
which  the  window  function  has  a  zero  in  its  Zak  transform,  this  being  one  of  the  singular 
situation  for  which  the  techniques  mentioned  above  were  devised.  The  philosophy  of  such 
methods  is  to  live  with  the  zero  rather  than  modify  it.  One  technique  is  to  represent  the  signal  as 
a  projection  into  the  space  of  all  signals  having  Zak  zeroes  where  the  window  function  has  them, 
and  minimizing  the  energy  loss  by  through  relative  positioning  of  the  signal  and  grid.  The  error 
involved  in  neglecting — in  Zak  space — a  coordinate  for  which  the  Zak  transform  of  the  window 
is  very  small  or  zero  is  considered,  and  the  minimization  of  this  error  over  shifts  of  the  signal  and 
window  is  shown  to  lead  to  smaller  than  errors  than  one  might  initially  expect.  A  specific 
example — white  gaussian  noise  input — is  carried  through  the  analysis.  Further  improvement  is 
possible  when  the  input  is  cyclostationary. 

We  also  indicate  how  a  slightly  modified  basis,  one  that  eliminates  a  number  of  the  Gabor 
translates  from  the  basis — as  many  as  there  are  Zak  zeroes  in  the  window — and  replaces  them 
with  other  functions,  provides  a  stable  basis  that  retains  most  of  the  Gabor  features. 

1.2.1.1J  Relations  between  Matrix  Stability  and  Signal  Duration 

In  section  2.1.3,  stability  of  the  matrix  method  of  Gabor  transform  calculation  is  studied  by 
looking  at  the  algorithm  that  produces  a  biorthogonal  function  from  its  window,  or  vice  versa. 
Our  preliminary  results  help  explain  the  mechanisms  controlling  the  blowup  of  the  transform 
method  for  certain  window  selections. 

1 .2.1.2  Algorithms 

1JL1.2.1  Oversampling  with  the  Gabor  Transform 

One  new  important  algorithm  has  been  implemented  in  GSPS  this  year,  which  is  the 
oversampled  Gabor  representation,  or  the  Weyl-Heisenberg  expansion.  In  this  theory  the  number 
of  expansion  coefficients  exceeds  the  number  of  data  samples,  which  leads  to  some  interesting 
tradeoffs.  For  the  higher  computation  required  to  obtain  the  coefficients,  one  gets  finer  tune- 
frequency  resolution — as  in  classical  Gabor,  one  can  partition  this  between  domains  at  will.  As 
could  be  expected,  this  induces  a  non  uniqueness  in  the  transform,  which  one  can  view  in  one  of 
two  equivalent  ways:  (1)  the  coefficients  associated  with  a  fixed  funcuon  are  not  unique,  or  (2) 
the  relationship  between  window  and  biorthogonal  is  no  longer  uniquely  invertable. 
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122  Software  Development 
1-2- 2-1  Second  Year  Effort 

The  accomplishments  of  the  second  year  effort,  during  which  the  bulk  of  the  GSPS  was 
developed,  are  reviewed  prior  to  discussing  the  third  year  amendments. 

i -2-2.2  Third  Year  Effort 

Among  the  principal  upgrades  of  the  GSPS  during  year  three  are:  a  map  view  of  the 
coefficients;  ability  to  turn  grid  lines  off  or  on  at  will;  improvement  of  the  command  line 
interface;  and  availability  of  intermediate  results  in  the  transforms.  In  accomplishing  the  latter, 
the  transforms  were  split  so  that,  for  example,  the  computation  can  be  halted  after  the  calculation 
of  the  7-aif  transform  without  continuing  through  the  Gabor  coefficient  caiculation.Also,  a  list  of 
warnings  is  provided  about  features  of  the  code  for  which  undebugged  traps  have  been  observed. 

1213  User's  Guide 

New  files  and  routines  are  listed,  and  as  an  example,  the  procedure  for  adding  a  new  signal  type 
is  explained.  Because  of  the  split-up  of  the  transform  calculations,  a  new  program  flow  is  in 
operation,  and  this  is  discussed  in  detail. 

1.23  Applications  Research 

i -2,3.1  Fault  Identification  in  Feedback  Control  Circuitry 

In  work  performed  during  the  first  two  years  of  the  contract,  AAEC  investigated  the  ability  to 
locate  faults  in  feedback  control  circuitry  using  data  sets  provided  by  a  DOD  customer.  In  this 
last  year  the  emphasis  was  on  the  test  of  morphological  filter  methods  for  this  problem,  under 
funds  nrirfryt  to  the  contract  through  a  MIPR.  A  significant  issue  in  the  evaluation  of  the  methods 
was  the  questionable  "ground  truth"  supplied  with  the  data.  Near  the  end  of  the  work  a  revision 
of  the  ground  truth  in  the  data  was  mad#.  using  a  neural  net  classifier,  and  subsequent  evaluations 
were  maria  against  this  standard.  Classification  using  the  methods  described  in  the  text  generally 
succeeded  in  achieving  a  high  correct  classification  probability. 

An  extension  of  the  Gabor  methods  used  on  this  problem  in  year  2  was  attempted,  but  the 
"ground  truth"  problems  associated  with  the  data  precluded  any  meaningful  effort,  since  the 
revised  ground  truth  was  not  available  in  time  for  the  needed  analysis. 


1-3 


Atlantic 


1 JL3 3  Automatic  Target  Recognition 

An  infrared  search  and  track  (IRST)  system  is  able  to  detect  targets  while  they  are  still 
unresolved  by  the  receive  optics.  Each  sufficiently  strong  target  generates  a  response  equal  to 
the  point  spread  function  (spatial  impulse  response)  of  the  optics.  The  premise  that  automatic 
detection  and  recognition  (ATR)  of  such  targets  in  a  background  of,  for  example,  clouds  could  be 
enhanced  by  Gabor-based  signal  processing  was  tested  under  contract. 

AAEC  employed  some  actual  cloud  background  data  upon  which  was  superimposed  both 
uniformly  and  nonuniformiy  distributed  targets  representing  point  spread  funcuons  of  either  a 
long  wave  or  medium  wave  IR  sensor.  Single  scan  lines  were  extracted  and  subjected  to  Gabor 
processing  characterized  by:  (1)  a  window  function  modeled  after  the  target  impulse  response, 
and  (2)  expansion  that  emphasized  time  (or  spatial)  resolution  at  the  expense  of  frequency 
resolution. 

The  experimental  results  showed  two  main  features.  The  first  is  that  if  one  performs  not  a 
single  Gabor  analysis,  but  instead  several  in  which  the  registration  of  signal  and  window  are 
changed  at  each  instance,  the  targets  essentially  all  show  up  in  one  of  the  cuts.  In  some  other 
cases  in  which  preprocessing  was  applied  to  the  signals,  the  targets  were  made  visible  in  a  single 
cut,  but  spread  across  the  frequency  bins. 

The  results  of  this  investigation  are  of  course  preliminary  and  incomplete,  but  they  do 
suggest  promise  for  the  Gabor  transform  as  a  tool  in  this  class  of  ATR  problems.  Applicauons 
involving  resolved  targets  were  not  tested. 

1.233  Minimum  Dimension  Gabor 

Under  a  prior  Phase  I  and  II  SBIR  for  ARP  A,  AAEC  has  performed  the  theoreucal  and  initial 
computational  exercises  in  minimization  of  the  dimension  of  a  representation  of  a  signal  set.  The 
efforts  carried  out  under  this  contact  further  generalize  this  by  incorporating  minimization  under 
linear  and  nonlinear  constraints  on  the  biorthogonal  function.  The  results  of  the  exercise  show 
this  process  to  be  tricky,  requiring  greater  effort  than  was  able  to  be  applied  here.  The 
computational  procedures  in  doing  this  have  by  now  been  pretty  well  wrung  out,  but  there 
remain  fundamental  questions  associated  with  the  best  types  of  constraints  to  use.  Also,  the 
ability  to  work  with  the  Zak  method,  which  we  believe  to  be  the  more  stable  of  those  available, 
has  vet  to  be  implemented. 
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CHAPTER  2 


THEORETICAL  DEVELOPMENTS 


2.1  THEORY 

2.1.1  Gabor  Sampdng  Theorem 

During  the  reporting  period  we  have  looked  at  the  application  of  sampling  theorems  within  the 
realm  of  Gabor  analysis.  This  activity  has  led  to  a  Gabor  sampling  theorem  that  casts  light  on 
the  relationship  between  degrees  of  freedom  in  the  Gabor  window  and  the  represented  signal. 
We  begin  with  a  naive  imitation  of  the  original  Shannon- KLotelnikov  theorem  derivation  and 
encounter  some  interesting  results  that  differ  from  the  classical  theory  due  to  the  supervening 
circumstance  that  there  are  two  functions  of  the  independent  variable  (taken  to  be  time  here) 
subjected  to  sampling.  This  imposes  some  stipulations  not  present  in  the  earlier  theory. 

To  begin,  we  address  the  use  of  Shannon-  Kotelnikov  sampling  expansions  for  functions 
expressed  in  a  Gabor  representation.  In  general,  a  function  represented  as 

(2.1.1) 

If Ul 

where  the  {wmjl(r)}  comprise  a  Gabor  basis  for  some  T  >  0, 

-  nDexpijlTtmt  /  T),  (2.1.2) 

does  not  admit  of  a  sampling  theorem  because  it  is  not  necessarily  bandlimited  Even  if  the 
window  function  w  is  itself  bandlimited.  (2.1.1)  can  have  arbitrarily  high  frequency  content  if  the 
Gabor  coefficients  do  not  vanish  as  m  increases.  But  with  a  handling  J  "nndow  and  Gabor 
coefficients  that  are  zero  beyond  a  maximum  frequency,  any  linear  combination  of  basis 
functions  will  be  bandlimited,  hence  subject  to  a  sampling  theorem.  We  develop  such  a  theorem 
and  use  it,  along  with  assumptions  of  time  limitation,  to  investigate  relationships  among  the 
time-bandwidth  products  of  the  signal  and  window  and  the  number  of  nonzero  time  and 
frequency  coefficients  in  the  expansion. 
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2.I.I.I.  The  Sampling  Theorem 

Let  w(t)  be  a  function  bandlimited  to  \f\ <  S*.  where  /  is  cyclic  frequency,  such  that  the  set  of 
translates  (2.1.2)  is  a  basis,  for  L2(R).  We  denote  the  set  of  such  handlimncd  functions  as  £{£*} 

and  say  that  we  £{5*}.  Correspondingly,  the  translate  e  B{BW  +  m/T}  Vne  Z.  For 

an  expression  of  the  form 

iW/2  « 

X  (2A3) 

mm  .Mt1nm~^m 


we  clearly  have  that  {h^}  eB{Bw  +  M / 27*}Vn  e  Z,m  £  M/1,  and,  by  Shannon- Kotelmkov 
[SHAN], 


wmji  (0  sinc|^2  7cB0  ^  t  -  j  j 

-  -  «rjexp(</2joiii  /  1BJ)  smcj^afl^f  -  ^  j  j, 


(2.1.4) 


where  B0  =  B^+M  /IT  and  we  use  sinc(«)  to  denote  sin  (•)/(•).  To  use  a  common  set  of 

samples  of  w  for  each  value  of  n  in  (2.1.4),  we  should  set  T  equal  to  a  multiple  of  1/2 BQ.  Le., 
T  =  Lq/  IBq.  Now  (2.1.4)  reads 


(2.1.5) 


Inserting  the  sampling  theorem  expansion  for  w  into  (2. 1 2)  yields 

s(t)~  X  I  amjt  £  ^:^^W/2^/A>)sinc[2xfi0f/-T|-)  .  (2.1.6) 

mm-M/lnm—.  i— —  V“o/  L\o/J 

Trivially  s  e B{B0 } ,  and  it  has  the  sampling  expansion 


(2.1.7) 


where 


'K)' 


(2.1.8) 


2.1.12.  Tbne-Bandwktth  Products 

No  restnction  on  the  rime  index  n  of  the  Gabor  coefficients  was  needed  to  get  the  sampling 
theorem.  But  now  let  us  apply  the  restriction  |nj<  N 1 2  and  consider  only  those  expansion  using 
a  finite  number  of  coefficients: 


M/2  Nil 

**)m  £  £<Vi.*Vn(0.  (2.1.9) 

m—M/2iv*—N/2 

Further  suppose  that  w  is  essentially  time-limited,  Le.,  the  samples  of  w  are  zero  for  |rj  £  Tw  /  2 . 
We  can  associate  with  w  a  time-bandwidth  product  (or  number  of  degrees  of  freedom) 

DW^2BWTW,  (2.1.10) 

that  represents  the  mininmum  number  of  samples  of  required  to  represent  w  faithfully  through 
the  sampling  theorem.  The  same  can  be  done  for  s,  and  it  is  easy  to  show  that 

Dsm2BsTs=(Tw  +  NT)(2Bw  +  M/T).  (2.1.11) 

Because  both  Z*>  (=  2 BoD  and  M  are  integers,  and  =  2BWT + M ,  the  quantity  2B#T  is  also  an 
integer,  which  we  denote  as  L,  the  length  of  the  Gabor  time  cell  in  units  of  Nyquist  samples  of  s. 
Then 


Dj  =*  (1  +  AT  /  L\2ByJw  +  NL).  (2.1.12) 

Assume  the  window  function  and  the  coefficient  bounds  M  and  N  are  given,  and  that  only 
the  grid  shape  (L)  may  be  varied.  To  varying  L  is  to  change  the  subspace  of  L2  ( R )  spanned  by 
the  basis  of  the  truncated  representation.  Although  we  would  like  our  basis  to  accommodate  as 
many  degrees  of  freedom  as  possible,  it  is  not  particularly  desirable  to  have  that  number  depend 
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on  the  grid  shape.  A  value  of  time-bandwidth  product  that  is  achievable  for  any  grid  choice 
would  be  a  conservative  guide  to  the  capability  of  the  truncated  representation. 

Minimizing  the  value  of  Ds  with  respect  to  U  however,  determines  a  value  of  signal  time- 

bandwidth  product  that  is  achievable  independent  of  the  choice  of  grid.  Signals  whose  time- 
bandwidth  products  are  less  than  this  minimum  will  be  called  uniformly  representable  with 
respect  to  the  window  and  the  coefficient  limits  { w,  M,  N) .  Ds  has  a  unique  minimum  at  L  a  L  , 


(2.1.13) 


(2.1.14) 


when  L*  >  1.  We  interpret  (2.1.14)  to  say  that  a  signal  is  uniformly  representable  by  the  basis 
{ w,  M,  N]  if  its  time-bandwidth  product  satisfies 

+JMN.  (2.1.15) 

If  (2.1.15)  does  not  hold  for  some  s,  then  either  the  time  extent  of  the  signal  will  exceed  the 
span  of  the  time  translates  of  w  or  else  s  will  have  frequency  content  that  is  not  captured  by  the 
frequency  translates,  at  least  for  some  choices  of  T  (or  I).  One  of  the  following  inequalities  will 
be  violated  for  some  choices  of  grid; 


Bs  >  BW+M/2T 

TS>TW  +  NT  (2.1.16) 

Although  no  theorem  guarantees  that  a  signal  uniformly  representable  with  respect  to  a  truncated 
basis  { w,  Af,  N)  will  automatically  have  a  ‘nice’  representation  in  the  sense  of  some  error 
measure,  for  example 


e  =  t(0-  X<W*m.n(')|  .  (2.1.17) 

S  112 

the  time-bandwidth  product  guidelines  should  provide  insight  into  the  match  between  a  signal 
and  its  expansion  over  a  truncated  basis.  Theorems  relating  to  this  will  be  sought  in  future  work. 
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We  close  by  observing  tint  (2.1.15)  can  be  solved  in  the  following  two  interesting  forms:  the 
first  states  a  requirement  on  the  window  time-bandwidth  product,  given  the  signal  and  the 
number  of  coefficients;  the  second  tells  how  many  coefficients  must  be  retained  to  expand  a 
given  function  with  a  particular  window: 

^2ByvTw  >^2BJS  -VStf 

4MN  t  <^22^57  -  tJIBwTw  •  (2.1.18) 

2.12  Accuracy  and  Stabilty  of  the  Expansion 

Accuracy  and  stability  of  the  Gabor  representation  are  longstanding  issues.  Over  the  course  of 
work  on  the  three-year  contract,  AAEC  has  made  considerable  progress  the  in  understanding  of 
these  matters,  most  of  which  has  been  discovered  piece  by  piece  and  not  documented  in  prior 
reports. 

Two  distinct  approaches  are  reported  here,  hr  the  first  type,  a  “bad”  window  is  tamed  by 
adjustments  to  either  the  window  or  the  point  grid.  Finding  a  transformation  close  to  the 
original,  but  with  much  better  conditioning,  is  the  goal  of  this  approach.  In  a  second  method,  we 
examine  the  implications  of  working  with  window  functions  that  lead  to  nonin vertible  transforms 
and  assessing  the  associate  loss  of  representation  fidelity.  The  latter  work  is  in  its  early  stages 
and  is  not  yet  supported  by  numerical  experiments,  whereas  the  former  category  has  been  well 
explored  numerically  by  AAEC. 

2.1JL1  Stabilization  by  Window  or  Grid  Adjustment 

In  [BALA],  presented  at  the  1992  IEEE-SP  International  Symposium  on  Time-Frequency  and 
Time-Scale  Analysis,  we  have  collected  descriptions  of  a  number  of  the  techniques  we  have 
successfully  employed  in  stabilizing  ill-conditioned  Gabor  expansions  and  presented  them  by 
computed  example.  Methods  applying  to  both  the  Zak  and  matrix  methods  are  included.  The 
following  summarizes  the  content  of  this  paper,  which  is  found  in  full  in  Appendix  B. 

The  paper  begins  by  summarizing  the  Zak  and  matrix  algorithms  used  in  AAEC’s  GSPS 
software,  highlighting  the  features  of  the  window  function  that  in  each  case  can  lead  to 
singularity  of  the  transform.  Section  2  addresses  experiments  performed  with  the  Zak  algorithm, 
where  we  present  a  condition  number  expression  for  the  mapping.  Explosion  of  this  number 
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results  from  zeros  of  the  Zaic  transform  of  the  window  function  lying  on  or  near  the  computation 
grid.  Using  a  window  of  the  forzn 

w(t)  *  _  -  (2.1.19) 

which  has  a  single  zero  at  a  point  whose  time  coordinate  is  adjustable  by  choice  of  a,  we  explore 
three  techniques  for  ameliorating  the  Zak  zero  problem:  (1)  subtle  change  of  the  window 
waveform  to  displace  the  Zak  zero;  (2)  tune-  or  frequency-translation  of  the  window;  and  (3) 
alteration  of  the  time-frequency  grid.  Examples  of  an  ill-conditioned  transform  before  and  after 
each  of  these  methods  are  shown,  and  in  each  case  considerable  stability  is  recovered. 

Section  3  deals  with  a  similar  treatment  of  the  matrix  method.  The  failure  mechanism  in 
this  case  is  the  occurrence  of  zeros  of  the  window  function  itself  within  the  interval  from  the 
zero-th  to  the  first  time  grid  point  The  window  in  (2.1.19),  having  a  zero  at  the  time  origin, 
makes  its  matrix  transform  nonin vertibie.  Left  shifting  the  window  by  just  one  point  is  enough 
to  stabilize  the  transform.  A  second  method  involves  adding  a  small  constant  value  to  the 
window  at  every  point  Numerical  results  of  both  techniques  are  presented. 

2.1^2  Working  with  Transforms  Having  Zak  Zeros 

The  following  thoughts  were  stimulated  by  reading  [TOLI],  in  which  it  is  remarked  that 
although  one  approach  to  deal  with  a  window  having  a  Zak  zero  is  to  use  it  only  for  expansion  of 
functions  that  have  a  corresponding  Zak  zero,  this  solution  may  be  overly  restrictive.  The 
reference  contains  no  evaluation  of  the  degree  of  restriction  imposed  by  this  approach,  and  it 
seemed  of  interest  to  consider  this  aspect  further. 

Assume  a  discrete  Gabor  expansion  using  a  window  having  exactly  one  Zak  zero  on  the 
grid,  £.g.,  the  gaussian,  for  which 

=  0  (2.1.20) 

In  a  discrete  time  expansion  there  is  just  one  bad  point  for  the  transform  as  executed  by  the  Zak 
method,  which  is  the  point  where  the  ''divide  by  zero'’  occurs.  Consider  the  set  of  all  functions  / 
in  l1  such  that 


Z/OS.ja-o. 


(2.1.21) 
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This  set  is  a  subspace  of  l2,  since  if  /  and  g  have  the  zero  property,  af  +  bg  has  it  also.  We 
would  like  to  consider  what  it  means  to  take  an  arbitrary  /  and  make  some  “minimal" 
modification  to  it  such  that  the  resulting  f  belongs  to  the  subspace  , %) »  0 . 

The  indiramd  procedure  is  of  course  a  projection  operator,  and  thus  there  is  a  least  squares 
solution  to  the  estimation.  That  is,  find  /  such  that:  (i)  and  (ii)  j/-/]2  is 

minimized.  Finding  /  could  be  accomplished  by  a  classical  least  squares  procedure,  but  there  is 
a  better  way.  One  should  compute  the  Zak  transform  of  /,  zero  the  appropriate  coordinate,  and 
re  transform  to  get  /. 


Let  h  represent  the  error  induced  by  the  projection. 


/  =  /  +  *. 


(2.1.22) 


and  use  as  the  error  metric.  Let  the  Gabor  expansions  under  consideration  have  an  N  x  M 
time  x  frequency  grid  size;  then  there  are  MN  points  in  any  associated  function  (the  data, 
window,  Gabor  coefficients  and  the  Zak  transforms).  In  all  the  following  we  assume  the 
gaussian  case  with  its  zero  at  the  center  of  the  unit  square.  It  is  easy  to  see  that 


f  0  ;  (/>.<?)*(%%) 

1 2(/)0W$>;  (*«>-(%  5$) 


(2.1.23) 


and  that 


fO; 

(-1)' 

N 


q*M/2 
q  =  M/  2' 


One  then  readily  finds  that 


(2.1.24) 


(2.1.25) 


(2.1.26) 


We  can  use  flAg2  /g/jpas  a  measure  of  energy  loss  in  the  projection,  and  estimate  its 
approximate  value  by  noting  that  the  loss  is  one  out  of  MN  coordinates:  since  the  Zak  transform 
is  unitary,  we  can  anticipate  an  energy  loss  on  the  order  of  UMN. 
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Consider  further  that  there  are  hi  unique  alignments  of  /  with  the  rjme  grid  in  regard  to  Zak 
transformation,  and  there  is  no  reason  to  imagine  an  a  priori  preference  for  any  one  of  them. 
Define  the  set  of  time  translates  of /as  {/^M>},  0  £  m  £  M  - 1,  where  we  identify /with/*°\  and 

/*m)  is  its  m-th  translate.  Then  in  attempting  to  minimise  the  projection  error,  we  are  free  to 
inspect  all  hi  translates  of  /  and  choose  the  one  such  that  jz(^l&))04>)^)j2 is  least.  The  error 
associated  with  each  translate  is 


f 

I 

jz(/<m,(%yM 

(2.1.27) 


and  the  minimum  is 


=  nun  4- 

(KmSiW-i 


(2.1.28) 


We  would  like  to  estimate  the  savings  in  choosing  to  discard  the  smallest  possible  amount 
of  energy.  To  do  so  requires  some  assumption  about  the  signal  being  represented;  let  us  as 
an  example  a  signal  /  consisting  of  statistically  independent,  real-valued,  zero  mean,  unit 
variance  gaussian  noise  variables.  For  any  fit  is  true  that 

ZfOi.%)’  £<-l)*/(*+K>.  (2-1.29) 

*-0 

and  therefore  the  new  gaussian  variable  in  (2. 1 .29)  satisfies 

Z^KJT)  =  0  (2.00) 

and 


| ZfM.yitf  -21 

km  0 


(2.1.31) 


We  would  like  to  compare  this  average  energy  in  the  Zak  transform  of  /  to  the  expected 
energy  in  the  Corresponding  Zak  of  signal  f(A),  for  which  the  energy  loss  is  minimum.  We  can 

do  this  through  the  following  theorem. 
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(2.132) 


K-  (M'^(  *7S- 

This  theorem  tells  us  that  there  is  an  expected  savings  proportional  to  the  square  root  of  the 
available  number  of  shifts  of  the  signal  relative  to  the  grid. 

Proof:  Let 

OZmZM-l  (2.1.33) 

and  define 

•S«min{|$J}.  (2.1.34) 

JVI 

We  want  to  find  the  probability  density  of  S  and  compute  its  second  moment 
Clearly 

Pr{S2X}  =  Pr{(k|>  X)n(|S,|>  X)— nflsw.,|>  X)}  =  [ft{(fS»|>  X)}]"  (ZJJJ) 

(2.1.36) 

(2.1.37) 

(2.138) 


But 

Pr{(N>  X)}  *  *{&  >  *M$>  <  -*)}  =  2o(i]. 

where  Q(*)  is  the  normal  probability  tail  integral 

Q(x)  *  J du  p(u)  m -1-  Jd«exp(-a2  / 2). 

X  X 

and  p(m)  is  the  gaussian  density.  We  observe  that 

4 -fioo =-/**)• 


The  probability  distribution  (integral  of  the  density)  of  s  is  then  given  by 


(2.139) 


and  its  density  is  found  by  diffexentiadon: 


*{{S  <  X)(  -  ;xiO 


(2.1.40) 


Then 


■  2  MN  j  dv  v2p(Vi2Q(v)f-1 . 

o 


(2.1.41) 


We  see  that  (2.1.41)  yields  a  value  of  N  for  the  case  M  »  1,  which  checks  with  (2.131). 
Although  we  cannot  analytically  cany  out  the  integration  in  (2.1.41)  to  get  a  closed  farm  exact 
value,  we  can  get  a  tight  upper  bound  by  using  the  familiar  inequality 

Q{v)  1  iexp(-v2 / 2);  v£0  (2.1.42) 


as  follows: 


c  dv_ 

J  llx 

0 

n)m&-  QED-  (ZL43) 
o 

Our  interpretation  of  this  result  is  as  follows.  If  we  arbitrarily  position  our  signal  with 
respect  to  the  grid,  we  can  expea  to  lose  the  energy  of  one  of  our  A fN  coordinates  If  we  replace 
the  numerator  and  denominator  in  (2. 1 .27)  with  their  expected  values,  we  find 


2  exp^-Afv2  /  2) 


Sz  » IMNjdv  v2p(vX2G(v)]*“‘  S  2 MN 
0 
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Choosing  the  signal  placement  wisely  lets  us  achieve 


(2.1.44) 


(2.1.45) 


Le.,  the  expected  fraction  of  energy  loss  is  smaller  than  the  fraction  of  coordinates  discarded. 
Thus  from  an  energy  viewpoint,  representation  by  a  discrete  Gabor  transform  using  a  window 
having  a  single  zero  may  be  quite  tolerable! 

There  is  another  perspective  to  put  on  this  sequence  of  observations.  We  have  selected  for 
our  signal  a  stationary  process,  and  have  applied  to  it  a  time-frequency  transform.  Conventional 
wisdom  tells  us  the  utility  of  time-frequency  techniques  is  that  they  capture  nonstationary 
behavior.  In  one  sense  then,  we  might  not  expect  any  benefit  from  the  analysis,  and  one  could 
argue  such  a  case,  saying  “our  window  selection  was  so  inept  that  the  basis  functions  were 
incomplete  and  we  lost  some  signal;  only  by  a  clever  trick  did  we  minimize  that  loss.” 

In  support  of  this  argument,  notice  that  the  expected  energy  loss  increases  as  W  increases, 
Le.,  as  the  time  resolution  of  the  transform  increases.  This  may  be  considered  to  be  the  penalty 
associated  with  going  from  Fourier-type  representations,  which  have  no  time  resolution,  to  a 
mixed  domain  picture  when  the  phenomenon  under  study  is  stationary.  Now  suppose  we  give 
our  signal  a  second-order  statistical  cyclostationarity  by  assuming  the  variance  to  be  a  periodic 
function  with  a  period  such  that  several  cycles  are  captured  within  the  MN  data  points.  In  the 
period  is  in  fact  M ,  we  would  anticipate  that  a  Gabor  representation  using  M  frequency  points 
and  N  m  MN  /  M  time  points  to  be  well  matched  to  the  signal.  If  the  variance  had  a  pattern  such 
as  that  shown  in  Fig.  2.1-1,  one  can  see  that  by  picking  the  grid  such  that  tire  time  coordinate  of 
the  zero  of  the  window  matches  the  minimum  variance  point,  one  would  in  fact  expea  to  lose 
little  and  do  better,  on  the  average,  than  the  1  /  4m  improvement  found  for  the  stationary  case. 
Here  we  are  exploiting  the  time-varying  behavior  of  the  Gabor  series  to  match  the  problem,  and 
are  finding  some  good  rewards. 
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Observing  ail  that  has  been  done  above,  one  can  conclude  that  it  is  not  actually  necessary  to 
yreiffrft  any  of  the  signal  energy,  if  we  will  liberalize  our  representation  slightly.  The  space  into 
which  we  have  projected  the  signal  has  an  orthogonal  complement  spanned  by.  in  this  case,  a 
cmgi-  basis  function.  If  we  add  this  function  to  the  (Af  N  -  1)  others  generated  by  the 
ume/firequency  translates  of  the  Gabor  window,  we  recover  a  complete  basis.  It  differs  from  a 
strict  Gabor  basis  in  that  one  member  fails  to  exhibit  the  translation  property.  This  gives  us 
another  way  to  see  the  penalty  for  using  the  window  with  a  Zak  zero;  it  causes  us  to  abandon  pan 
of  the  nice  structure  we  valued  in  going  to  the  Gabor  series  in  the  first  place.  However,  the 

deviation  can  be  considered  minimal.  In  fact,  the  added  basis  function  depends  only  on  one  of 
«irh  M  points  in  the  data;  this  function  is  just  g**h/ Zf  (%,%),  where  h  is  given  by  (2.1.22); 


g 


0; 

-  (-1 Y 

.  N 


q*M/ 2 
\q  =  MH 


(2.1.46) 


Observe  g  is  zero  except  for  one  point  in  M. 

The  above  viewpoint  says  that  we  need  not  abandon  a  window  function  that  resembles  our 
data  quite  well  simply  because  of  its  Zak  zeros.  The  options  are  to  remain  within  the  strict  Gabor 
structure  and  sacrifice  a  small  amount  of  signal  energy,  or  to  enlarge  the  representation  and 
capture  everything.  One  of  these  approaches  may  well  suit  a  wide  variety  of  problems. 
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2.13  Relations  between  Matrix  Stability  and  Signal  Duration 

The  Matrix  Gabor  representation  as  described  in  [BALA2]  has  some  stability  problems 
associated  with  it  that  ate  mentioned  in  the  original  paper,  such  as  nomnverdbility  of  the  window 
if  it  posesses  an  analytical  (or  computational)  zero  in  the  first  Gabor  time  slice,  as  well  as  other 
problems  that  axe  mentioned  in  the  previous  section  and  were  investigated  in  fBALAl],  a  copy  of 
which  is  included  in  Appendix  A.  We  will  now  describe  some  preliminary  results  which  lead  us 
to  believe  that  one  may  have  some  control  over  the  stability  of  the  method  by  controling  the 
signal  length  or  the  number  of  Gabor  time  points,  depending  on  the  nature  of  the  signal  that  is  to 
be  analyzed.  For  windows  whose  values  in  the  first  Gabor  time  slice  are  not  zero,  but  have  a  big 
dynamic  range  and  with  the  smallest  value  being  much  less  than  one,  the  method  is  also  unstable 
as  can  be  seen  from  the  following.  Let  the  linear  system  to  be  solved  be  defined  as 

(WE)x=*b  (2.1.47) 

where  W  is  the  MN  x  MN  (block  )  matrix  containing  the  window  values,  £  is  the  MN  x  MN. 
(diagonal)  Fourier  rotation  matrix,  x  is  the  MN  X  1  vector  of  (ordered)  unknown  coefficients,  and 
b  is  the  MN  x  1  vector  of  signal  data  points,  M  is  the  number  of  Gabor  frequency  points  and  N  is 
the  number  of  Gabor  time  points.  We  can  now  write  the  vector  x  of  coefficients  as 

x  =  E'Bb  ,  B-W“  (2.1.48) 


Let  us  now  write  W  and  B  in  their  block  form  in  order  to  see  their  representation  more  clearly 
and  to  see  how  B  will  depend  on  the  structure  of  W.  Note  that  because  of  the  ordering  each  of 
the  blocks  of  B  and  W  will  be  M  x  M,  and  there  will  be  N  x  N  of  them  in  each  matrix  (for  a  full 
derivation  see  (BALA2]). 


■  *0 

*0 

w  = 

*7 

. 

*2 

V. 

"V.1 

... 

w. 

(2.1.49) 
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b,  =.-wjl(w,w;1) 

B,  - -Wj'fWjW;1  +  W,B,  +  W,B,) 


aijo) 


(2.1J1) 


B,  -  -WJl(W,W^  +  WwB,+-+W1B ,_t);  V/  e  [ZN - 1] 

The  extent  of  our  analysis  to  date  has  been  to  look  at  the  simplest  nontrivial  case  which 
occurs  when  the  window  lasts  only  over  the  first  two  Gabor  time  slices.  This  being  the  case,  the 
zeroth  and  the  first  W  are  nonzero,  and  the  above  formula  for  the  hiarthogonai  matrix  becomes 

Bo— W? 

B1  a-W^fWjWJ1) 

B,  =  -W^l(W1B1)  =  W^1(W1W^1)2  (2U2) 

B3  =  -  W^1  ( WjB.  )  =  -  W^1  (Wj  W^1  )3 


B,  =-W^l(W1B/_1)  =  (-l)iW^1(W1W^1y;  V/e[2.1V-l] 

from  which  we  notice  two  important  things.  The  first  is  that  even  though  the  window  only  has 
support  over  two  time  slices,  the  biorthogonal  has  support  over  all  the  time  slices,  therefore,  the 
biorthogonal  matrix  will  be  full  More  apropos  to  this  discussion,  we  notice  that  each  successive 
block  of  the  biorthogonal  grows  by  a  power  of  the  inverse  of  the  window  entries  in  the  zeroth 
time  slice  multiplied  by  the  corresponding  window  entries  in  the  first  time  slice.  This  is  due  to 
the  diagonal  nature  of  each  of  the  blocks  of  W  and  B.  This  implies  that  the  stability  of  the 
method  depends  on  containing  that  growth  rate  to  reasonable  bounds. 
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The  above  discussion  implies  that  one  now  has  two  possible  ways  of  condoling  the  stability 
of  the  method  depending  on  the  nature  of  the  signal.  If  the  signal  of  interest  is  such  that  it  can 
tolerate  being  split  up  into  pieces,  then  one  can  choose  the  desired  number  of  frequencies  M  and 
then  make  the  length  of  the  signal  be  such  that  the  n»mher  of  time  points  N  keeps  the  power 

Bv_,  *(-l)*“‘W?(W1W?)*"‘.  (2.1.53) 

well  behaved  and  bounded.  One  can  then  apply  the  transform  to  each  of  the  pieces  separately.  If 
on  the  other  hand,  the  signal  in  question  does  not  tolerate  being  split  up,  but  we  have  some 
freedom  in  choosing  the  number  of  frequency  points,  ie,  the  window  is  nonzero  over  the 
increased  length  of  the  Gabor  time  slice,  then  one  can  increase  the  number  of  frequency  points 
thereby  decreasing  the  number  of  time  points  (for  a  fixed  signal  length)  and  keeping  the  above 
expression  well  behaved. 

The  computation  of  the  biorthogonal  matrix  entries  gets  more  involved  as  the  number  of 
Gabor  time  points  that  support  the  window  inmeases,  and  this  analysis  has  not  yet  been 
performed  for  the  general  case.  When  the  analysis  is  complete,  we  will  submit  a  paper  for 
publication  on  the  obtained  results.  As  a  closing  note  we  remark  that  in  GSPS.  all  the 
aforementioned  analysis  can  be  done  a-priori  for  windows  of  all  lengths  of  interest,  and  the 
information  will  be  readily  available  when  the  signal  is  chosen,  since  the  analysis  is  completely 
independent  of  the  signal. 

22  ALGORITHMS 

22.1  Oversampling  with  the  Gabor  Transform 

In  order  to  increase  the  capability  of  the  GSPS  workstation  to  allow  us  to  analyze  different 
types  of  signals,  we  decided  to  include  the  capability  to  perform  oversampling.  A  brief 
description  of  the  algorithm  is  given  here,  and  for  a  full  detailed  explanation  the  reader  is 
referred  to  [WEXL].  The  formulation  of  the  oversampling  algorithm  is  as  follows:  Let  P  be  the 
total  number  of  points  in  the  signal  data  set,  M  and  N  be  the  number  of  desired  frequency  and 
time  points  respectively,  and  M'  and  AT  be  two  auxiliary  positive  integers  satisfying  the  constraint 
that 


P  =  MN*  +  M'N 


(2.2.1) 


Under  these  conditions,  the  Gabor  transform  can  be  written  as 


(2-L2) 


M-M-l 


sik)»  XZ<~  h(k-mtr)cxQ(2ixnk/ N) 

r-\ 

*  ^s(k)b‘(k - mAOexp(2iJDik  /  AO 


(2*2*3) 


together  with  the  condition  that 


+  mAOA  (*)exp(2iJDiJ: /  AA)»(P/  MN)SmS, . 


(2JL4) 


At  this  stage,  the  last  condition  can  be  rewritten  in  matrix  form,  in  order  to  solve  for  the 
biorthogonai  function  b  in  a  similar  way  as  the  matrix  formulation  of  the  Gabor  transform.  The 
M *  x  M  block  matrix  of  window  values,  each  square  block  being  AT  x  N' .  can  be  written  as 

(WE)b*r  (2J2S) 


(WE)b-r 

*0 

Wvr 

^iur-OAr 

W 

^■urHU-Dtr 

"(tf-iW**’ 

W 

"(w-o tr*iir 

or 

-m/h-oat 

with  each  Wj  given  by 


^h*or-D 


it 


(22.1) 


for  a  grand  total  of  M'  WxP  entries,  premultiplied  by  the  Af ’  xM’.  block  matrix  Eo  of  Fourier 
rotations  of  size  AT  given  by 
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E< 


X 


X 


X 


(2-2.8) 


and  the  P  x  1  unknown  biortho  gonal  and  the  M'  W  x  1  right  hand  side  vectors  b  and  r  are 
respectively. 


b*[l*0)  1*1)  1*2)  b(P-2)  i*P- l)]r 

r  *  [P  /  MN  0  0-0  0]r 


(2.2.9) 


The  biortho gonal  vector  can  be  now  calculated  after  choosing  a  way  to  'invert1  the  non-square 
matrix  associated  with  the  window  function.  We  chose  to  implement  a  generalized  inverse  or 
energy  minimization  method  in  the  GSPS  workstation  with  some  slight  computational 
modifications  from  the  referenced  paper.  We  noticed  that  since  the  inverse  is  given  by 

b  «  (EW)r((EWXEW)r)Hr 

3Wr(WWr)-f  ,  (2JL10) 

?  =  EHr 


the  last  computation  can  be  performed  by  hand  and  it  results  in  a  reasonable  computational 
savings.  The  new  right  band  side  vector  becomes 


r»[r0  0  0  -  0  Of 

r0=[P/MN  P/MN  P/MN  -  P/MN]t' 


(2J1.11) 


with  each  one  of  the  blocks  being  of  size  N'  x  1  .  The  biorthogonal  is.  at  this  stage,  fully 
available  and  the  Gabor  coefficient  map  can  be  computed  with  the  aid  of  (2.2.3).  Reconstruction 
of  the  signal  (with  or  without  postprocessing)  can  now  be  also  ratod?"**  with  the  aid  of  (2,2,2) 
A  better  description  on  how  to  perform  these  operations  within  GSPS  and  a  description  of  which 
parameters  can  be  chosen  by  the  user  will  be  described  in  Section  3.3. 
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CHAPTER  3 

SOFTWARE  DEVELOPMENT 


3.1  SUMMARY  OF  SECOND  YEAR  EFFORT 

Before  discussing  the  software  effort  expanded  this  third  and  last  year  of  the  contract,  we 
will  briefly  summarize  the  work  that  was  performed  during  the  second  year  under  contract. 

3.1.1  Phase  n  Code 

During  the  second  year  of  effort,  numerous  additions  and  enhancements  were  incorporated 
into  the  Gabor  Signal  Processing  System  software,  which  from  here  on  we  will  call  GSPS. 
Before  moving  on  to  summarize  the  Phase  II  code,  we  will  begin  by  describing  the  philosophical 
change  that  took  place  towards  the  end  of  the  first  year  of  the  contract.  As  described  in  the  First 
Annual  Technical  Report,  the  initial  Gabor  transform  software  prototype  was  implemented  on  an 
IBM  compatible  PC  and  it  was  command  line  driven,  ie,  the  user  was  prompted  to  enter  a 
□umber  corresponding  to  a  signal,  window,  etc.,  the  Gabor  processing  was  performed,  and  the 
results  were  then  written  to  files.  Unfortunately,  in  order  to  change  parameters,  one  had  to  either 
restart  the  code  or  hardwire  different  values  into  the  program  (such  as  total  number  of  points)  and 
recompile.  The  cumbersomeness  of  doing  this,  together  with  memory  limitations  and  speed 
considerations,  forced  us  to  migrate  to  a  different  (better)  platform  and  to  create  a  user  friendly 
graphical  user  interface  that  could  accomodate  all  the  power  and  flexibility  that  we  wanted  to 
incorporate  into  the  GSPS.  Bearing  that  in  mind,  at  the  beginning  of  year  two  we  selected  a 
SUN  Microsystems  SP ARCStation  1  +®  as  the  computational  plaform,  C  as  the  scientific 
programming  language,  and  XI 1  together  with  the  XView  toolkit  as  the  graphical  user  interface 
programming  language  to  be  used  to  develop  the  GSPS  software  that  is  currently  available,  and 
which  we  have  been  using  in  house'  for  a  period  of  well  over  one  year. 

3.1.1.1  Functionality 

The  organization  of  the  program  has  become  more  clear  to  the  user  by  the  development  of  a 
graphical  user  interface  that  is  layed  out  in  such  a  way  that  the  user  can  follow  a  sequence  of 
major  categories,  represented  on  the  screen  as  butttons,  and  progresivciy  move  down  the  logical 
sequence  to  the  desired  function  to  be  applied.  This  is  accomplished  in  the  following  manner 
the  major  categories  represented  by  buttons  are: 


Signals 

Windows 

Noise 

Method 

Clipping 

Reconstruction 

Iterate 

Utility 

Coefficients 

Options 

Output 

Exit  GSPS  ? 

Each  of  the  buttons  is  the  owner  of  a  menu  which  contains  menu  items,  which  in  turn  may 
contain  (nested)  submenus.  The  last  entry  of  a  path  through  a  menu  is  either  a  command  which  is 
executed  immediately  or.  when  appropriate,  a  command  which  contains  a  dialog  box  in  which 
the  user  can  either  enter,  change,  or  choose  selected  default  parameters.  These  parameters  are 
then  updated  in  the  main  program  and  the  command  is  executed.  If  the  action  that  is  selected  by 
the  user  results  in  a  graphical  representation  of  data,  a  window  is  automatically  opened  and  the 
Hata  is  displayed.  Since  the  real  estate  on  a  computer  screen  is  limited,  all  of  the  windows  that 
are  opened  by  GSPS  are  multifunctional,  le,  the  same  graphical  window  that  is  used  to  display 
the  analysis  window  is  used  to  display  the  function  that  is  biorthogonal  to  said  analysis  window. 

Keeping  that  in  mind,  we  can  now  illustrate  the  high  level  execution  of  the  program  by 
m^nnc  of  an  example,  while  a  more  detailed  example  including  the  features  incorporated  in  the 
third  year  will  be  included  in  the  User  Manual  (Section  3.3).  The  user  starts  execution  of  GSPS 
by  opening  their  favorite  XI 1 -based  window  manager  on  the  screen,  and  typing  gsps  from  one 
of  the  text  I/O  windows,  or  alternatively,  opening  a  file  manager  and  double  clicking  on  the  gsps 
application.  This  results  in  the  display  of  a  Control  panel  containing  the  (  deactivated  except  for 
Exit  GSPS)  buttons,  and  a  Disclaimer  panel  containing  the  proprietary  information,  disclaimer, 
software  version,  and  two  buttons  marked  "Continue"  and  "Quit".  Selecting  the  Quit  button  exits 
GSPS  with  no  action  being  taken,  while  selecting  the  Continue  button  erases  the  Disclaimer 
panel  and  activates  the  buttons  in  the  Control  panel.  We  are  now  ready  to  start  processing 
signals.  Selecting  the  "Signals"  button  with  the  Right  Mouse  Button  (  RMB  )displays  a  menu 
containing  different  signal  choices  right  under  the  Signal  button.  Selection  of  a  signal  choice 
with  the  Left  Mouse  Button  (LMB)  opens  up  a  dialog  box  which  allows  the  user  to  either  choose 
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default  values  for  the  signal  by  selecting  tbe  “Defaults"  button  in  the  dialog  box,  or  enter  the  user 
defineabie  parameters  manually.  Tbe  action  can  now  be  either  aborted  by  selecting  tbe  "Cancel" 
button,  or  accepted  by  selecting  the  "Ole”  button.  Accepting  the  choice  results  in  the  appearance 
of  a  display  window  which  contains  the  selected  signal.  One  can  now  proceed  to  the  selection  of 
the  analysis  window  by  selecting  the  "Window"  button  with  the  RMB.  This  results  in  the  display 
of  a  menu  that  is  functionally  identical  to  the  "Signals"  menu  and,  in  fact,  contains  the  same 
entries.  Selecting  the  "Ok”  button  from  the  Dialog  box  results  in  the  appearance  of  a  second 
display  window  that  contains  the  graphical  representation  of  the  analysis  window.  Once  the 
signal  and  window  have  been  chosen,  the  user  can  proceed  to  the  choice  of  the  method  that  will 
be  used  to  do  the  processing  by  chosing  a  method  from  the  menu  associated  with  the  "Method" 
button.  This  can  be  accomplished  by  selecting  the  "Method”  button  with  the  RMB.  and  choosing 
a  method  with  the  LMB  results  in  a  third  display  window  being  open,  and  it  contains  the 
appropriate  coefficients  corresponding  to  the  selected  method.  It  is  noteworthy  to  mention  that 
the  Coefficients  window  automatically  displays  both  2-D  and  3-D  coefficient  sets  depending  on 
the  representation  that  is  warranted  by  the  given  method.  At  this  stage  the  user  can  choose  to 
repeat  the  portion  of  the  experiment  described  so  far  by  choosing  different  signals,  windows, 
methods,  or  any  combination  of  the  above,  or  proceed  to  the  reconstruction  of  the  signal.  If  one 
opts  for  the  reconstruction  route,  there  is  an  option  to  first  perform  simple  thresholding  on  the 
coefficient  set  by  selecting  the  "Clipping"  button  with  the  RMB,  entering  the  desired 
thresholding  levels,  and  removing  the  coefficients  below,  above,  or  between  user  defined  upper 
and  lower  bounds.  If  it  turns  out  that  it  is  not  desireable  to  perform  any  thresholding,  one  can  go 
directly  to  one  of  the  reconstruction  routines  which  behave  analogously  to  the  method  options  in 
the  "Method"  button.  Tbe  reconstructed  signal  is  now  displayed  in  the  Signal  window  and.  if  the 
reconstruction  method  is  the  same  as  the  transform  method,  the  coefficient  set  has  not  been 
clipped,  and  the  method  was  stable,  it  will  be  the  same  as  the  original  signal  except  for  roundoff 
error.  The  1*2  distance  between  the  original  signal  and  the  reconstructed  signal  is  reported  next 
to  the  origin  in  the  Signal  display  window 

3.1J  Capabilities 

During  the  second  year  of  effort,  numerous  additions  and  enhancements  were  made  to  the 
tool  that  was  available  at  the  end  of  the  first  year,  the  most  noteworthy  of  which  being  that  the 
user  could  change  initialization  parameters  like  signal  duration  and  time-frequency  splitting  of 
the  grid,  as  well  as  other  parameters,  without  having  to  recompile  the  program.  In  fact,  the  only 
time  that  recompilation  is  needed  is  when  a  new  feature  is  added  to  one  of  the  menus  as  a  menu 
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item.  Alter  the  addition  is  debugged,  ail  the  parameters  that  can  be  con  trolled  by  the  user  (if  any) 
will  appear  in  a  dialog  box  when  the  selection  is  made. 

Another  major  improvement  was  the  feel  of  the  tool.  The  user  coultToow  input  parameters 
of  interest  in  the  dialog  boxes,  and  choosing  the  button  labeled  'Ok1  made  the  program  accept  all 
the  parameters  at  once  and  continue  with  the  execution.  If  the  wrong  parameter  was  input,  the 
user  could  move  the  pointer  over  to  the  appropriate  field  and  change  it  before  selection  of  the 
'Ok'  button  if  the  mistake  was  noticed  in  time  or,  in  the  case  that  it  was  not.  the  same  (or  a 
different)  menu  item  could  be  chosen,  the  correct  parameter  input  and  the  'Ok'  buton  selected.  In 
the  previous  version  of  the  software,  clerical  errors  would  result  in  the  user  having  to  restart  the 
program  fro  step  one.  The  menu  approach  also  allowed  the  user  to  pick  up  at  any  point  on  the 
list  of  buttons  which  logically  preceded  the  last  step  chosen,  and  after  an  initial  run  through  the 
four  logically  sequential  steps,  ie,  after  going  through  the  process  of  choosing  the  signal, 
window,  method  and  reconstruction,  one  could  reprocess  from  any  of  the  intermediary  steps.  We 
will  illustrate  by  saying  that  after  the  user  had  chosen  the  signal,  window,  and  method,  the  signal, 
window  or  method  could  be  chosen  again,  as  well  as  the  associated  utilities,  without  having  to 
either  complete  the  sequence  by  choosing  the  reconstruction  option  or  starting  out  with  a  new 
copy  of  the  same  (or  different)  signal.  The  logical  flow  of  the  program  is  illustrated  in  Fig.  3.1*1 
below. 


Fig.  3.1*1.  Logical  flow  of  the  GSPS  package  at  the  end  of  the  second  year 


Note  that  even  though  the  user  could  choose  a  different  window  or  method  after  having 
gone  through  the  sequence  of  events  once,  the  grid  definition  was  at  the  signal  level,  therefore,  if 
one  wanted  to  redefine  the  grid  shape  and  size,  one  would  have  to  go  back  to  the  signal  choice 
and  redefine  the  signal  in  such  a  way  that  it  reflected  the  desired  configuration.  We  will  see  in 
the  next  section  that  this  limitation  no  longer  exists  due  to  the  ability  of  the  window  to  redefine 
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grids,  and  the  added  capability  to  redefine  grid  configuration  from  a  menu  item's  dialog  box 
under  the  Utility  button. 

During  The  second  year  of  effort  we  also  included  the  method  resulting  from  the  matrix 
reformulation  of  the  transform  [BALA]  and  numerous  utilities  to  preprocess  the  signal,  such  as 
zooming,  shifting,  saving  to  memory,  etc..  We  also  reported  various  values  of  interest  to  the 
appropriate  window  after  certain  actions  were  taken,  such  as  the  length  of  the  maximum 
coefficient  and  the  dimension  of  the  signal  set  [ORR?].  For  a  foil  description  of  the 
enhancements  introduced  during  the  second  year  the  reader  is  referred  to  [AAE1]. 

33  THIRD  YEAR  EFFORT 

3.2.1  New  Facilities 

As  reported  in  the  previous  technical  report  at  the  end  of  year  one  it  was  decided  to  improve 
our  computational  facilities  in  order  to  be  able  to  process  larger  data  sets,  and  to  speed  up  the 
processing  time.  As  it  turns  out  the  upgrade  that  was  made  to  the  system  at  that  time  has 
resulted  in  more  than  adequate  computational  facilities  to  perform  the  tasks  asociated  with  this 
project  and.  therefore,  no  new  facilities  were  needed  during  the  third  year  of  effort.  For  a  foil  list 
of  facilities,  the  reader  is  referred  to  [AAEC]  and  [AAE1]. 

333  Phase  HI  Code 

Phase  m  software  development  has  centered  around  broadening  the  applicability  of  GSPS 
to  a  wider  range  of  problems  by  improving  its  ability  to  interoperate  with  other  analysis 
packages,  and  through  the  addition  of  new  diagnostic  displays.  In  addition,  we  have 
implemented  two  new  versions  of  the  Gabor  transform  with  slightly  different  capabilities  than 
those  present  in  the  Phase  II  code.  To  facilitate  to  the  addition  of  these  features  to  GSPS,  several 
key  portions  of  the  Phase  n  code  have  been  restructured  for  added  flexibility,  notably  in  the  file 
I/O  facilities  and  in  the  handling  of  signal  buffers. 

333.1  New  Display  Components 
map  view 

A  “map”  mode  has  been  added  to  the  coefficient  display  to  supplement  the  spike  plot 
display  from  Phase  II.  In  this  mode,  the  coefficient  matrix  is  displayed  as  a  grid  of  rectangular 
color  patches  whose  colors  correspond  to  the  magnitude  of  the  coefficient  at  that  coordinate, 
relative  to  the  maximum  coordinate  magnitude.  Coefficient  values  below  a  compiled-in  epsilon 
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value  are  suppressed,  and  are  displayed  in  the  background  color.  A  color  key  is  provided  to  the 
lower  left  of  the  grid.  Selection  of  the  spike  plot  and  map  modes  is  accomplished  using  the 
“map”  button  located  in  the  control  area  above  the  color  key. 

grid  button 

A  button  labeled  “grid”  now  appears  in  the  control  areas  of  the  signal,  window,  and 
coefficient  displays.  In  the  signal  and  window  displays,  this  option  controls  the  display  of  the 
hash  marks  which  indicate  the  endpoints  of  the  Gabor  tunesiices.  In  the  coefficient  display,  this 
option  works  in  conjunction  with  the  map  button  to  control  the  display  of  reference  lines  in  the 
spike  plot  mode,  or  to  mark  the  boundaries  of  each  color  patch  in  the  map  mode. 

The  RMS  duration  and  RMS  bandwidth  are  now  computed  for  both  the  signal  and  window, 
and  are  displayed  in  the  control  areas  of  their  respective  windows. 

A  number  of  new  signal  options  ate  available: 

(interaction  changes:  new  signal  and  window  types) 

Get  Fourier 
null  signal 

(interaction  changes:  reading  data  files) 
getfileinfo,  matlab 

(interaction  changes:  new  transforms) 
matnx2,  oversampling 

(interaction  changes:  piecewise  forward  and  reverse  transforms) 

the“compute”  routines 

bioithogonais  and  windows:  matrix  and  zak 

forward  and  reverse  zak 

(interaction  changes: ) 

(changes  to  utilities) 

One  .significant  departure  from  the  Phase  II  code  is  that  the  buffers  containing  the 
intermediate  components  of  the  Gabor  transforms  (i.e.  the  biorthogonai  and  Zak  buffers)  may 
now  be  loaded  from  a  file  and  explicitly  manipulated.  This  introduces  a  synchronization 
problem  as  these  buffers  were  implicitly  linked  to  other  buffers  when  their  contents  were 
automatically  gemrated.  A  set  of  global  flags  have  been  added  to  indicate  the  validity  of  these 
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buffers,  and  should  be  checked  before  the  corresponding  buffers  are  to  be  used.  Likewise,  these 
flags  should  be  set  appropriately  whenever  any  of  the  buffers  is  modified. 

3JL3  Warnings  and  Suggestions 

The  following  is  a  list  of  known  caveais  for  potential  users. 

After  the  tool  has  been  running  for  a  long  time,  it  may  crash. 

Since  you  can't  'edit  data  in  the  ‘window’  window,  any  Gabor  window  preprocessing  can  be 
done  in  the  signal  window,  copied  to  memory,  and  then  read  in  the  window  window. 

The  number  of  points  in  a  signal  or  window  should  be  an  integer  power  of  2. 

Although  the  coefficient  display  is  cleared  whenever  it  becomes  marked  as  invalid,  a  refresh 
call  to  that  window  may  cause  the  previous  contents  to  be  redisplayed. 

A  number  of  conventions  must  be  observed  when  Matlab  data  files  are  to  be  imported  into 
GSPS.  Several  variable  names  are  reserved  for  use  by  GSPS  [tmpM  tmpN  tmpRate  tmpType 
tmpMethod  Comment].  Any  stored  data  elements  which  do  not  correspond  to  one  of  these 
names  will  be  interpreted  as  the  data  with  which  the  buffer  will  be  filled.  As  a  result,  you  should 
not  save  any  other  elements  in  this  file  except  for  the  signal  itself,  which  may  be  saved  with  an 
arbitrary  name.  If  the  file  importing  routine  encounters  an  unknown  element  following  the 
signal,  it  will  attempt  to  load  that  element's  contents  instead. 

GSPS  was  developed  in  the  Open  Windows™  2  environment,  and  is  designed  primarily  for 
use  in  a  color  environment  Although  the  application  will  run  successfully  with  a  monochrome 
display,  various  text  fields  will  appear  crowded  as  the  monochrome  XView™  2.x  libary  uses  a 
different  font  entirely.  GSPS  will  compile  and  run  with  XView  3.x  libraries  as  well,  but  will 
exhibit  this  behavior  on  both  color  and  monochrome  displays. 

3.2.4  Scripts  ' 

Although  the  Phase  Q  version  of  GSPS  included  a  simple  command  line  interface,  it  has 
since  proven  inadequate  to  the  expanded  capabilities  of  the  Phase  HI  code.  As  GSPS  has  become 
more  general  in  its  purpose,  so  has  the  need  for  a  more  generic  command  line  interface  grown. 
The  new  interface  syntax  is  designed  to  accommodate  additional  transforms  and  also  makes 
provisions  for  transform-specific  parameters. 
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gsps  -s  <fnam*>  -w  <fname>  -xf2lmioizJ  -o  <fnam*> 

gsps  - s  <fname>  -w  <jhame>  -m  <M>  -ft  <N>  -xf2\miolz}  ~o  <fname> 

when: 

•s  <fname>  is  a  path  to  the  file  containing  the  signal 
-w  <fname>  is  a  path  to  the  file  containing  the  window 
x  specifies  the  transform  to  be  used: 

-x2  indicates  the  Manix2  method 
xm  indicates  the  Matrix  method 
•xo  indicates  the  oversampling  method 
-xz  indicates  the  Zak  method 

-o  <fname>  contains  the  name  of  the  file  to  which  the  coefficients  should  be  written 
33  USERS  GUIDE 
3.3.1  Procedure  Descriptions 
The  following  new  files  and  routines  are  present: 

new  files:  get_four_coef.c  getfileinfo.c  loarimat.c  null.sig.c  over_gabor.c  savemaix  tbprodx 
uttlmatx 

new  routines:  graphicsx  initx  ioadfile.c  mainx  savefile.c  xforms.c 

In  the  event  that  a  signal,  window,  or  coefficient  matrix  of  exceptionally  small  magnitude 
must  be  displayed,  the  plotting  routines  now  impose  a  non-zero  scale  value... 

Example:  how  to  add  a  new  signal  type 

Although  new  signals  are  typically  generated  externally  and  imported  through  the  file  I/O 
interface,  it  is  frequently  desirable  to  hardcode  new  signal  types  which  are  to  be  used  frequently 
with  minor  variations,  even  though  these  signals  will  only  be  accessible  from  the  graphical  user 
interface.  This  procedure  typically  requires  four  steps:  coding  the  signal  generator,  coding  the 
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user  interface  components,  adding  calls  to  the  main  initialization  and  menu  generation  routines, 
and  adding  build  dependencies  to  the  compilation  control  configuration  file.  For  clarity,  it  is 
generally  desirable  to  name  the  source  code  files  after  the  signal  being  generated;  for  instance,  a 
triangular  pulse  generator  might  have  the  files  tri_pulse .  c  and  optionally,  tr impulse .  h 
associated  with  it.  Because  GSPS  lacks  a  generic  mechanism  for  manipulating  signal  generation 
parameters  outside  of  the  GUI.  these  variables  should  remain  opaque  to  the  remainder  of  the 
program  and  may  be  defined  as  static  types  local  to  the  signal’s  source  module. 

The  signal  generation  routine  is  often  adapted  from  another  piece  of  standalone  software,  so 
typically  it  is  passed  a  desired  number  of  points,  a  sample  rate,  and  whether  a  real  or  complex¬ 
valued  signal  is  desired,  and  typically  it  returns  a  buffer  of  the  appropriate  size  and  type.  This 
may  then  be  called  from  a  callback  routine  which  is  shared  by  the  signal  and  window  menu 
items,  which  is  then  responsible  for  copying  the  generated  signal  into  the  correct  buffer. 

The  GUI  supporting  code  consists  of  an  initialization  routine  which  creates  the  GUI  objects, 
and  a  set  of  callback  routines  which  those  objects  will  dispatch  as  various  buttons  and  widgets 
are  manipulated.  The  construction  of  these  callbacks  is  beyond  the  scope  of  this  document,  but 
the  file  tri_pulse .  c  is  recommended  as  a  template. 

Modifications  should  be  made  to  the  file  main .  c  in  four  places.  First,  external  references 
to  the  routines  should  be  added  near  the  top  of  the  file.  Next,  menu  items  should  be  created 
under  both  the  Signal  and  Window  menus.  Finally,  a  call  to  the  initialization  routine  should  be 
added  near  the  end  of  the  file,  but  before  the  main  loop  is  invoked. 

3-3.2  Program  Flow 

During  the  third  year  effort,  much  functionality  was  added  to  the  GSPS  software,  but  the 
logical  flow  of  the  program  has  changed  very  little  since  the  end  of  the  previous  year.  The  main 
logical  change  has  been  due  to  the  inclusion  of  a  menu  item  called  grid,  info  to  the  Utility 
button.  This  menu  item  allows  the  user  to  redistribute  the  time  and  frequency  points  on  the  grid 
at  any  time  after  the  signal  buffer  has  been  written  to  at  least  once  has  been  added  to  the  utility 
button,  therefore,  the  user  no  longer  needs  to  choose  a  new  (or  the  same)  signal  to  run  different 
grid  configurations.  The  user  still  needs  to  generate  a  new  signal  buffer  if  the  signal  duration  is 
to  be  changed.  One  other  change  worth  mentioning  is  that  under  the  Utility  button,  a  Compute 
menu  item  has  been  added  which  allows  the  user  to  compute  intermediate  steps  in  the  Gabor 
transform  computation,  for  instance,  the  Zak  transforms  and  the  window  biorthogonal.  This 
means  that  the  user  no  longer  has  to  perform  the  full  transform  computation  if  the  only  thing  of 
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interest  at  any  particular  tune  is,  say,  the  biorthogonals  to  the  window  that  result  from  a  set  of 
gndspadngs.  For  big  signals,  this  can  prove  to  be  a  considerable  tune  savings 
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CHAPTER  4 

APPLICATIONS  RESEARCH 

4.1  FAULT  IDENTIFICATION  IN  FEEDBACK  CONTROL  CIRCUITRY 

This  section  describes  the  experiments  performed  to  investigate  the  application  of  the  Gabor 
representation  for  solving  problems  of  emitter  identification  in  communications.  In 
experiments,  test  data  was  used  to  investigate  the  characterization  of  signal  features,  specifically 
transient  signals,  with  the  Gabor  transform.  These  test  data  features  were  used  to  characterize 
signals  as  belonging  to  specific  classes  of  emitters  as  an  exercise  in  determining  utility  of  Gabor- 
based  methods  for  this  problem 

The  problem  investigated  in  these  experiments  was  the  mapping  of  signals  of  interest  to 
transmitting  and  receiving  equipment  identity.  The  utility  of  using  Gabor-based  rerhmqiM»«  was 
studied  and  tested  with  a  set  of  signal  data  from  several  transmitter-receiver  pairs.  pair  was 
considered  to  generate  a  class  of  signals.  The  Gabor  techniques  were  applied  to  the  dam  to 
characterize  each  signal  as  a  member  of  one  of  these  classes  of  signals.  A  dimension  measure  of 
the  Gabor  coefficient  set,  developed  by  AAEC  under  a  separate  DARPA  contract  [AAEC2],  was 
the  main  discriminant  tested  to  separate  the  signals  into  classes. 

4.L1  Summary  of  Previous  Work  (Greenbelt  Facility) 

4.L1.1  Data  Sets 

Two  sets  of  signal  data  were  used  during  this  effort  The  first  database  consisted  of  186  files 
containing  approximately  65  thousand  time  samples  each,  sampled  from  analog  signals  at  a  20 
Khz  rate.  Of  the  186  files,  only  47  files  contained  signals  of  interest  Two  types  of  transient 
features  were  found  in  these  47  files,  and  56  occurrences  of  these  features  were  chosen  for  the 
experiments.  The  two  types  of  transient  features  were  designated  the  small  and  normal  size 
signals.  The  length  of  these  features  of  interest  ranged  from  256  to  2048  time  samples. 

A  second  database  consisted  of  1944  files,  containing  signal  data  from  many  different 
transmitter-receiver  pain.  The  length  of  these  files  ranged  from  approximately  1000  to  16500 
time  samples.  In  this  data  set,  a  signal  of  interest  contained  about  200  time  samples. 
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4.1.122  Experiments 

The  technical  approach  used  in  this  investigation  was  to  visually  examine  the  original  signal  data 
and  look  for  features  in  the  signals  which  appeared  to  be  common  within  each  class  of  signals. 
Once  these  types  of  features  were  ide?  tified,  a  portion  of  the  data  containing  the  feature  was 
extracted  to  be  used  as  both  a  signal  and  as  a  Gabor  window  function.  Data-denved  window 
functions  were  then  applied  to  the  other  signals  using  the  Gabor  transform.  The  resultant 
dimension  of  the  set  of  Gabor  coefficients  was  minimiMri  through  a  time  alignment  procedure, 
and  this  minimum  dimension  was  used  to  discriminate  between  classes.  Additional 
discriminants  were  used  as  appropriate.  These  included  the  maximum  amplitude  of  the  Gabor 
coefficients,  and  the  reconstruction  error  due  tc  "clipping"  of  the  Gabor  coefficients,  which  is 
defined  as  deletion  of  all  coefficients  of  magnitude  smaller  than  a  percent  of  the  largest 
magnitude  coefficient.  These  additional  measures  improved  the  classification  process. 

Three  types  of  experiments  were  performed:  1)  Normal  vs  small  signal  charactenzaiion.  2) 
Small  signal  discrimination,  and  3)  Signal  classification  with  averaged  windows.  Since  signals 
of  interest  in  the  data  come  in  two  classes,  the  first  task  was  the  discrimination  between  these, 
addressed  by  the  first  experiment.  Given  success  at  this  first  stage,  discriminauon  within  a  class 
was  the  remaining  key  factor.  Experiments  2  and  3  were  devoted  to  discrimination  among  small 
signals.  Discrimination  among  the  normal  size  signals  was  not  emphasized. 

4.1.123  Results  and  Observations 

Even  though  only  a  finite  number  of  test  signals  and  six  window  functions  were  used  to  analyze 
and  characterize  the  signals  of  interest,  use  of  the  minimum  Gabor  dimension  value  and  the 
maximum  Gabor  amplitude  showed  promising  results.  The  use  of  averaged  windows  was  found 
to  be  a  useful  additional  method  of  discrimination.  For  a  complete  description  of  these 
experiments,  refer  to  the  Second  Annual  Technical  Report  (AAEC2). 

4.122  Summary  of  Previous  Work  (Waltham  Facility) 

During  the  previous  year,  a  MIPR  subprogram  was  conducted  under  the  Gabor  program  to 
address  the  problem  of  identifying  specific  emitters  from  a  particular  class  of  interest  to  the 
MIPR  sponsor.  A  number  of  very  large  (megabyte-size)  signals  were  provided  by  the  MIPR 
sponsor.  Each  typically  contained  several  transient  signals  of  roughly  a  thousand  samples  each. 
The  sponsor  provided  the  locations  of  a  few  of  these  transients.  The  remainder  (which  comprised 
the  majority  of  the  available  transients)  were  extracted  by  hand  by  plotting  the  million-sample 
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signals  as  1000  x  1000  "images"  and  manually  marking  interesting  segments  with  an  interactive, 
cursor-based  utility  which  we  developed  for  this  purpose.  Later  in  the  program,  a  morphological 
algorithm  for  automatically  detecting  and  locating  transients  was  developed  and  partially  tested. 
On  the  MIPR  sponsor's  request,  work  on  automatic  detection  and  extraction  was  terminated  so 
that  the  efforts  could  be  concentrated  on  emitter  classification 

Morphological  filtering  comprised  the  core  of  the  classification  techniques  developed  and 
tested,  particularly  for  signal  conditioning  and  feature  extraction.  The  development  of  reliable, 
quickly  computable,  morphology-based  discriminants  and  means  for  utilizing  them  was 
continued  into  the  third  year  of  the  Gabor  program,  and  is  discussed  in  the  next  section. 

4.1.3  Year  3  Effort  (Waltham  Facility) 

A  second  MIPR  was  added  to  the  Gabor  program  to  continue  work  on  the  specific  emitter 
identification  problem.  The  primary  technical  objective  was  to  demonstrate  the  effectiveness  and 
efficiency  of  morphological  and  other  modem  signal  processing  techniques  for  specific  emitter 
classification  using  transient  signals  provided  by  sponsor.  The  four  key  steps  performed  towards 
this  end  were  the  preprocessing  of  selected  transients,  feature  extraction,  classifier  design,  and 
performance  evaluation.  Morphological  techniques  were  applied  during  the  first  two  steps,  and 
were  found  to  have  the  most  benefit  for  feature  extraction,  providing  a  compact,  easily  computed 
representation  of  the  transients. 

Two  data  sets  were  provided  by  the  MIPR  sponsor  during  Year  3,  hereafter  referenced  as 
the  old  and  new  data  sets,  respectively.  The  old  data  set  comprised  23  moderately  long  signals, 
typically  containing  1000  samples  each.  Since  our  mam  thrust  was  classification,  i.e.,  not 
transient  detection  and  location,  we  manually  extracted  a  single  segment,  typically  100  sample 
long,  from. each  signal.  Although  the  Year  3  "old  data  set"  signals  were  much  smaller  than  those 
from  the  previous  year,  they  still  required  manual  transient  extraction  since  they  also  contained 
other  signal  components  which  the  sponsor  specifically  instructed  us  to  ignore.  The  signatures 
taken  from  the  old  data  set  were  fairly  clean.  The  ground  truth  provided  with  them  divided  them 
into  six  classes,  and  appeared  reasonable  based  on  both  visual  and  computer-aided  analysis. 

The  new  data  set,  also  provided  by  the  (MIPR)  sponsor,  was  much  larger,  consisting  of  1944 
moderately  long  signals  of  roughly  1000  samples  each.  180  segments,  typically  100  samples 
long,  were  manually  extracted  from  the  collection.  In  contrast  with  the  old  signatures,  many  of 
the  new  ones  had  a  very  noise-like  appearance,  probably  due  to  distortion,  dispersion,  widely 
varying  band-limiting,  and  other  effects  induced  by  the  transmission  channel  through  which  the 
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signals  were  received.  Furthermore,  it  was  difficult  to  relate  the  ground  truth,  which  segregated 
the  signatures  into  1 1  classes,  to  the  observed  waveform  characteristics. 

Similar  experiments  were  conducted  on  both  of  the  aforementioned  data  sets.  A  high  degree 
of  success  was  achieved  with  the  old  data:  all  but  one  of  the  23  signals  were  classified  correctly. 
The  lessons  learned  in  both  defining  the  architecture  of  the  various  classifier  components  and 
fine-tuning  the  associated  algorithmic  parameters  were  then  applied  to  the  new  data  Although  a 
number  of  different  feature-space  representations  were  tested  and  a  moderate  degree  of  success 
was  obtained,  it  became  clear  towards  the  end  of  the  second  MIPR  research  program  that 
information  on  the  underlying  physics  of  the  various  emitting  devices  (unavailable  during  either 
MIPR  program)  would  be  required  in  the  future  to  yield  reliable  classifier  performance. 
Nonetheless,  it  was  demonstrated  that  even  for  a  data  set  such  as  the  new  one,  morphological 
measures  provide  performance  comparable  to  or  better  than  conventional  techniques,  and  they  do 
so  with  less  demanding  computational  requirements. 

The  old  data  set  consisted  of  23  signals  collected  from  six  different  emitters.  The  signals 
from  one  of  the  six  classes  are  shown  in  Figure  4.1-1.  It  is  important  to  note  that  the  transients, 
although  similar  in  appearance,  are  never  in  any  way  time-aligned  or  "registered”.  The  classifier 
processing  to  be  described  is  completely  insensitive  to  the  transient  starting  points  within  the 
extracted  segments.  (Similar  insensitivity  can  be  achieved  by  conventional  means  such  as  ratring 
the  magnitude  of  the  transient's  Fourier  transform,  but  at  a  much  greater  computational  cost). 

A  high-level  block  diagram  of  the  classifier  system  is  shown  in  Figure  4.1-2.  As  mentioned 
above,  morphological  and  other  advanced  signal  processing  techniques  were  applied  in  the  first 
two  blocks.  Variants  of  an  efficient  traditional  classifier  were  used  for  the  last  block.  For  signal 
conditioning,  two  fundamental  morphological  filtering  operators,  opening  and  closing,  were  used 
to  remove  spikes  and  other  undesirable  temporal  characteristics.  Opening  consists  of  a 
morphological  erosion  followed  by  a  dilation.  Closing  consists  of  the  concatenation  of  these  two 
operators  in  the  reverse  order. 

A  two-dimensional  example  of  erosion  is  shown  in  Figure  4.1-3.  Like  all  basic 
morphological  operators,  erosion  is  based  on  a  structuring  kernel  which  is  analogous  to-a  finite 
impulse  response  (FIR)  filter.  The  kernel  is  “slid”  across  an  input  signal  or  image  just  as  is  done 
with  an  FIR  filter  during  convolution,  except  that  each  output  sample  is  obtained  from  the 
minimum  or  maximum  (for  erosion  and  dilation,  respectively)  of  the  input  samples  under  the 
sliding  kernel,  rather  than  the  weighted  sum  of  those  samples  las  in  convolution). 
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Fig.  4.1-2.  CJMtefler  system  block  diagram. 
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switches  shown  in  optimal  positions 

Fig.  4.1-4.  Signal  conditioning  experiments. 
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Many  combinations  of  morphological  and  more  traditional  signal  processing  techniques 
were  applied  to  the~raw  signals  prior  to  feature  extraction,  as  suggested  in  Figure  4.1-4. 
References  to  number  of  points  indicate  FIR  filter  or  morphological  kernel  sizes.  As  mentioned 
earlier,  the  segment  extraction  was  performed  manually.  The  switches  in  the  block  diagram  are 
shown  in  the  positions  which  ultimately  yielded  the  best  classification  performance.  They  show 
that  the  best  signal  conditioning  was  achieved  using  a  6-point  moving  average  filter. 


Fig.  4.1-5*.  Effect  of  signal  conditioning;  before  conditioning. 
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Fif.  4.1-Sc.  Effect  of  signal  conditioning;  as  input  to  feature  extractor. 
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The  effect  of  the  preferred  signal  conditioning  method  is  shown  in  Figure  4.1-5.  The  first 
two  plots  show  the  extracted  raw  signals  for  a  single  class  before  and  after  conditioning.  Note 
that  the  various  signals  in  theses  two  plots  have  been  time-registered  and  have  had  their  means 
removed  for  visual  comparison  only.  The  third  plot  shows  the  conditioned  signals  as  actually 
input  to  the  feature  extractors. 


Given  that  linear  filtering  was  judged  to  be  best  for  signal  conditioning,  the  dominant  role  of 
morphological  processing  was  in  feature  extraction.  For  this  purpose,  variants  of  a  higher-level 
morphological  product  known  as  pattern  spectra  were  used.  Like  openings,  closings,  and  other 
morphological  operators,  pattern  spectra  can  be  applied  to  data  with  any  number  of  dimensions 
(in  our  case,  just  one).  They  are  based  on  sequences  of  increasingly  larger  structuring  kernels 
which  are  matched  to  expected  waveform  features.  These  kernels  are  successively  applied  to  a 
given  signal  and  the  change  in  area  under  the  signal  stored  as  a  pattern  spectrum  "bin"  value  (see 
Figure  4.1-6).  By  convention,  the  change  in  area  due  to  larger  and  larger  openings  are  plotted  as 
positive  values  on  the  right  side  of  a  two-sided  pattern  spectrum.  Negative  values  are  plotted  on 
the  left  for  the  corresponding  closings.  (Figure  4.1-6  only  shows  the  openings  spectrum.). 


Filtering  Kernels 


Fig.  4.1-6.  Calculation  of  pattern  spectrum. 


Pattern  spectra  have  a  number  of  attributes  which  are  relevant  to  signal  source  identification. 
They  are  sensitive  to  small  signal  substructure  such  as  ringing  and  overshoot,  and  are 
unhampered  by  repetitive  patterns  in  the  data.  They  are  also  insensitive  to  macroscopic  attributes 
such  as  temporal  registration  and  overall  signal  amplitude.  Like  most  morphology-based 
products,  pattern  spectra  can  be  easily  and  rapidly  computed  in  real-time  using  special-purpose 
hardware  based  on  parallel  architectures,  and  do  not  suffer  from  the  dynamic  range  growth 
associated  with  most  conventional  techniques. 
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Figure  4.1-7  summarizes  the  various  feature  extraction  methods  which  were  applied  to  the 
"old"  data  set  (A  few  additional  techniques  were  applied  to  the  new  data  set.  as  discussed  later.) 
A  new  type  of  pattern  spectrum  was  developed  during  the  M1PR  program  and  found  to  be  the 
best  method,  among  all  those  tried,  for  characterizing  emitter  signatures.  It  employed 
logarithmically  spaced  pattern  spectrum  (LSPS)  bins  for  which  kernel  size  was  increased 
exponentially  (rather  than  linearly,  as  done  in  the  fine-grain  (FGPS)  and  coarse-grain  (CGPS) 
pattern  spectra).  The  log-spaced  pattern  spectra  for  the  signals  from  Figure  4.1-1  are  shown  in 
Figure  4.1-8. 

Once  candidate  features  were  extracted  from  all  test  signals,  they  were  fed  into  a  traditional 
maximum-likelihood  classifier.  Such  classifiers  assume  multi-variate  Gaussian  feature 
distributions  characterized  by  hyper-elliptical  contours  of  constant  probability  for  each  class  i,  as 
shown  in  Figure  4.1-9.  In  general,  each  estimate  of  the  class  S\omega_iS  to  which  a  given 
feature  vector  most  likely  corresponds  is  based  on  the  value  of  i  which  maximizes  the  log- 
likelihood  function 


(4.1.1) 
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For  simplicity,  we  assumed  that  the  a  priori  class  probabilities  Sp(\omegaJ)S  and  the 
covariance  matrix  determinants  St\Sigma_ilS  were  equal,  so  that  maximizing  the  log-likelihood 
.  function  from  above  reduces  to  mitiimmng  the  aonnalized-axes  distance  function 

(4.1.2) 

The  minimum  distance  (M3ND1S)  classifier  was  trained  in  three  different  ways,  as  depicted  in 
Figure  4.1-10. 


r 


switches  shown  in  optimal  positions 

Fig.  4.1-10.  Emitter  dassificattea  experiments. 


Mean  and  median  feature  vectors  were  found  to  be  equally  effective  in  characterizing 
specific  emitter  classes.  Means  are  somewhat  more  attractive  from  a  computational  standpoint 
since  they  do  not  require  sorting.  Simply  taking  the  first  feature  vector  from  each  class's  training 
set.  although  computationally  inexpensive,  was  (not  surprisingly)  less  effective.  The  results  of 
apply  the  mean-  (or  median-)  and  first-vector  training  methods  to  eight  different  types  of  feature 
vectors  are  summarized  in  Tables  4.1-  1(a)  and  4.1 -1(b).  Only  14  features  were  necessary  to 
achieve  the  best  performance  of  96%  correct  (highlighted  by  a  box  in  the  table)  using  the  mean- 
based  classifier.  Hie  other  boxes  in  the  table  highlight  signals  which  were  misclassified. 
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Table  4.1-1  •  NONDIS  classifier.  remits  for  old  data  set  (a)  Using  mean  or  median  feature  vectors  for  each 

dass.  (b)  Using  first  feature  vector  for  each  dass. 
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Versus  signal  from  true  class  (mmmlmmmh,  from  class  0): 


(•) 


-4  0  -20  0  20  <0  40 

sample  \pmmm  w  omo  loyauwon} 


(b) 


FT*.  4.1-11. 


i  typical  ripal  from  true  class,  (b)  Vc 


4-14 


Pina 1 TTfcAPPUCATlON  OF  THE  GABOR  REPRESENTATION  TO  MILITARY  PROBLEMS 


The  misclassified  4%  of  the  old  data  set  comprised  a  single  signal  from  among  all  23. 
Figure  4.1*11  provides  a  visual  comparison  of  this  signal  ("mmm3mmme",  from  Class  0)  with 
typical  signals  from  its  true  class  ("mmmlmmmh",  also  from  Class  O)  and  the  estimated  class 
("bs_lbs_a",  from  Class  A).  The  "time -constants'1  and  other  waveform  structure  associated  with 
"mmmlmmme"  are  clearly  much  more  similar  to  those  in  "bs_ibs_a"  than  those  in 
"mmmlmmmh". 


Table  4.1-  2.  MINDIS  rlantflur  reea Ha  for  extended  old  data  set. 


Morphological  processing  is  attractive  for  the  emitter  classification  problem  not  only  due  to 
its  capability  for  shape-oriented  signal  discrimination,  but  for  its  ease  of  implementation  and 
computational  efficiency  as  well.  To  evaluate  the  computational  complexity  of  our  techniques 
and  compare  them  with  more  traditional  methods,  let  us  define  the  following: 


•  the  signal  length,  L  =  71  to  501  (typically  approx.  100) 

•  the  number  of  bit  per  sample,  B  =  16 

•  the  moving  average  window  length,  Af  s  6 

•  the  maximum  log-spaced  pattern  spectrum  kernel  size.  K  =  64 

•  the  number  of  classes,  Q  -  6 

Let  us  assume  that  an  addition,  subtraction,  minimum,  or  maximum  computation  has  a  relative 
computational  cost  of  one  operation,  while  a  multiplication  requires  B  operations  (essentially  B 
additions).  The  unnormalized  short-term  moving  average  used  for  signal  conditioning  requires 

C„  =  (Af  -  l)L  adds  »  ML  ops  C4  i 
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For  each  kernel  size,  k,  the  log-spaced  pattern  spectrum  needs 

cn  k  a (k- 1  )L  min«  +(k  -1)L  maxs  +  2L  subs  +2(L-1)  adds  »2JkL  ops.  (4.1.4) 

Thus,  the  total  computatioa  for  the  log-spaced  pattern  spectrum  is 

C.  =  Yc^  -  V  2(2'  L)  *  2(2‘*,,(rHI  -  1)L  -  4 KL  ops 

*-i*«-*  «o  (4.1J) 


The  MINDIS  classifier  requires 

Cd  =2(logj(K)  +  l)Q  oiuls  +  2((log,(K)+l)-l)Q  adds  +(fi-D  mins 

*  2QBlog2(K)  ops  (4.1.6) 

The  total  computational  cost  is  thus  approximately 

C»(M  +  41T)£  +  2QBlog2(K)  ops/signal  (4.1.7) 

For  the  values  enumerated  earlier,  the  total  computation  is  dominated  by  feature  extraction,  i.e., 

C  -  4  XL  ops  *  27  kops  /  signal  .  (4.1.8) 

For  comparison,  the  computation  associated  with  Fast  Fourier  Transform  (FFT)  based  feature 
extraction  is 

Cm  •  Llogj(L)  cplx  -  mul  -  adds  =  4Llog,(L)  mul  -  adds 

-4flLlogj(L)  ops  (4.1.9) 


For  the  values  above, 

crrr  =  43  kops  /  signal  .  (4.1.10) 

which  compares  favorably  with  the  27  kops/signal  required  for  the  morphological  approach. 

In  summary,  the  performance  of  the  linear  filtering  signal  conditioning,  morphological 
feature  extraction,  and  minimum-distance  classification  worked  well  on  the  old  data  set.  The 
preferred  system  correctly  classified  22  (96%)  of  the  23  signals  in  the  old  data  set.  Comparable 
results  (97%  correct)  were  achieved  after  also  including  12  new  sponsor-provided  signals.  All 
processing  elements,  comprising  unnormalized  short-term  moving  averages  for  signal 
conditioning,  log-spaced  pattern  spectra  for  feature  extraction,  and  the  MINDIS  classifier  based 
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on  class-wise  vector  means,  were  computationally  inexpensive.  The  morphological  transform 
typically  reduced  signature  data  by  90%  and  obviated  the  need  for  time  registration. 

r  vea  the  success  above,  similar  techniques  were  applied  to  the  new  darn  set  Unfortunately, 
the  class  "truth"  provided  with  the  new  data  appeared  questionable  based  on  inspection  of  the 
plotted  waveforms:  signal  structure  varied  widely  between  signals  in  each  class  and,  more 
importantly,  there  were  many  suspicious  similarities  between  signals  in  different  classes  We 
therefore  decided  to  reassign  the  data  to  a  new  set  of  classes  using  sophisticated,  well-known 
classifier  tools  employing  neural  network  analysis,  and  then  apply  our  classifier  system  to  the 
data  with  the  new,  neural  net  truth.  In  addition  to  assigning  the  signals  to  new  classes — 
discarding  sponsor-defined  class  assignments — the  neural  network  program  was  intended  to 
determine  the  true  number  of  classes  represented  by  the  new  dam  set.  Once  the  classes  were 
reassigned,  the  preferred  preprocessing  and  classification  algorithms  which  were  successful  on 
old  data  were  applied  to  the  new  data,  and  alternate  feature  extraction  methods  were  developed 
and  tested. 

Neural  network  analyses  are  categorized  by  two  basic  paradigms,  supervised  and 
unsupervised.  Supervised  neural  nets  learn  and  generalize  existing  classes  and  are  most  suitable 
when  the  class  exemplars  are  well  known.  In  sharp  contrast,  unsupervised  neural  nets  are  used  to 
discover  classes  in  unclassified  data  and  thus  are  suitable  when  the  underlying  classes  are 
undetermined  or  vague,  as  was  the  case  for  the  new  data  set. 

The  neural  network  program  which  we  selected  for  use  as  our  "higher  authority"  for  class 
reassignment  was  SOM-PAK  (the  "Self-Organizing  Map  Program  Package"),  one  of  the  most 
readily  available  unsupervised  clustering  networks.  It  was  built  by  the  team  which  developed  the 
theory  [KOHO]  at  the  Helsinki  University  of  Technology  Laboratory  of  Computer  and 
Information  Science.  Its  one  disadvantage  is  that  it  does  not  automatically  estimate  the  number 
of  classes  in  the  data  and  thus  requires  a  hypothesized  number  of  classes  as  an  input. 

The  process  by  which  we  reassigned  classes  using  SOM-PAK  and  evaluated  our  own 
classifiers  is  summarized  in  Fig.  4.1-12.  As  evident  in  the  figure,  only  log-spaced  pattern  spectra 
were  input  the  neural  net  classifier,  whereas  all  types  of  feature  vectors  were  used  lor  testing  our 
classifier. 
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The  preliminary  experiment  conducted  using  SOM-PAK  was  to  search  for  any  degree  of 
correlation  between  the  sponsor-defined  and  neural  net  class  assignments  for  the  new  data  set.  It 
was  hypothesized  that  the  number  of  classes  was  20,  a  number  somewhat  larger  than  the  number 
of  sponsor-defined  classes  (11),  with  the  assumption  that  excess  classes  would  be  sparsely 
populated.  Table  4.1-3  consists  of  tallies  of  sponsor-defined  class  members  in  each  neural-net- 
defined  class.  There  is  no  apparent  correspondence  between  the  two  sets  of  class  assignments. 
Had  there  been  any,  there  would  have  been  only  one  relatively  large  value  in  each  row  (the  first 
row  contains  an  1 1.  a  6.  and  a  5;  the  second  contains  a  13  and  a  9;  etc.). 
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Table  4.1-3.  Tallies  of  sponsor-defined  dan  members  in  each  oenni-net-defined  class. 


The  next  preliminary  experiment  was  to  estimate  the  number  of  classes  spanning  the  new 
data  set  manually,  with  the  aid  of  SOM-PAK.  Table  4.1-4  illustrates  the  distribution  of  the 
signal  population  over  neural  net  classes  for  varying  hypothesized  numbers  of  classes.  For 
comparison,  boxes  are  shown  around  values  from  Table  4.1-3.  In  viewing  the  number  of 
members  per  class,  we  are  looking  for  a  cutoff  (a  row  number)  below  which  the  remaining 
classes  are  lightly  populated.  There  appear  to  be  two  dominant  (heavily  populated)  classes,  and 
about  five  classes  total  plus  individual,  random  (unassociable)  cases. 

In  the  absence  of  clear  information  on  the  number  of  underlying  classes  in  the  new  data  set, 
we  decided  to  evaluate  classifier  performance  for  various  numbers  of  classes,  as  shown  in  Table 
4.1-5.  The  values  shown  are  percent  correct  using  various  types  of  normalized  pattern  spectra 
(log-spaced,  etc.)  with  training  based  on  mean  feature  vector  codewords.  Classes  were  reassigned 
by  running  the  neural  net  on  the  LSPS  data  only  once  per  hypothesized  number  of  classes.  The 
best  performance  (76%  correct)  was  obtained  by  concatenating  the  LSPS  and  FGPS  feature 
vectors,  for  a  total  of  42  features. 
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Tabic  4.1-6.  Exploration  <rf  alternate  feature  extraction  methods. 


Alternate  feature  extraction  methods  woe  also  investigated.  The  original  (i.e.,  sponsor- 
provided)  truth  was  used  for  this  purpose  to  provide  a  fair  comparison  between  the  LSPS,  which 
had  also  been  used  for  class  reassignment,  and  other  types  of  features.  Table  6  compares  the  size 
and  effectiveness  of  the  various  feature  vectors.  All  pattern  spectra  vectors  were  unnormaiized. 
As  with  the  old  data  set,  the  best  performance  (38%  correct  for  the  new  data  set)  was  achieved 
with  the  pruned  LSPS  comprising  only  14  features.  Other,  much  more  computationally  intensive 
signature  measures  such  as  FFTs  and  parametric  (LPQ  modeling  did  not  improve  performance. 

4.1.3  Year  3  Effort  (Greenbelt  Facility) 

4.13.1  Experiments 

In  year  3,  the  second  database  of  emitter  data  was  processed  with  the  Gabor  software.  This  data 
consisted  of  1944  files  which  contained  several  instances  of  data  per  file.  Many  GSPS  runs  were 
made  using  signal  data  and  other  functions  as  windows.  In  these  cases,  the  minimum  dimension 
resulting  from  the  Gabor  processing  was  higher  than  in  previous  runs,  because  the  second 
database  contained  much  more  complicated  signals.  The  identification  of  which  portion  of  each 
file  was  the  signal  of  interest  was  difficult,  if  not  impossible  without  knowledge  of  the  source  of 
the  data. 
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After  many  attempts  at  classifying  elements  of  the  second  database,  it  was  determined  that  the 
second  set  of  data  was  unreliable,  and  the  "ground  truth"  could  not  be  known  due  to  errors  made 
by  another  contractor  when  generating  the  data. 

4-1-5  Future  Work 

It  is  clear  that  future  work  on  the  specific  emitter  identification  problem  will  be  most  effective  if 
the  physics  of  the  emitting  devices  are  analyzed  and  quantified.  This  was  not  possible  during  the 
programs  discussed  above  due  to  the  limited  availability  of  supplemental  data.  It  will  probably  be 
necessary  to  characterize  the  transmission  channel  through  which  the  signals  were  obtained,  and 
perhaps  the  data  acquisition  process,  as  well. 

Armed  with  all  of  this  information,  it  should  be  possible  to  implement  various  rapid  signal 
classification  processing  techniques.  Work  performed  using  the  Gabor  expansion  signal 
dimension  was  promising,  but  was  not  further  testable  under  this  year's  effort.  Extensions  might 
combine  Gabor  methods  with  those  based  on  pattern  spectrum  and  other  morphological 
discriminants.  The  excellent  results  (96%  correct)  obtained  for  the  old  data  set  of  23  signals 
suggest  that  the  log-spaced  pattern  spectra  will  provide  an  effective,  efficient  method  for  feature 
extraction,  in  terms  of  both  classification  reliability  and  feature  vector  compactness.  Other 
morphology-based  features  should  also  be  investigated  since  they  can  be  tuned  to  signal 
"appearance",  Le.,  they  can  discriminate  various  shapes  in  transients  in  a  manner  similar  to  what 
human  analysts  do.  Low  cost,  real-time  implementations  are  also  possible,  and  should  be 
developed  and  analyzed  with  respect  to  size,  cost,  and  power  requirements. 
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422  Problem  Definition 

42.1.1  Statement  and  Importance  of  the  Problem 

The  detection  of  small  targets  in  clutter  is  a  problem  of  critical  importance  to  wide-area,  long- 
range.  InfraRed  Search  and  Track  (IRST)  surveillance  applications.  In  these  applications,  the 
targets  (aircraft  and  missiles)  are  typically  unresolved  and  therefore  appear  in  sky,  sea  and  terrain 
backgrounds  in  only  a  few  resolution  cells.  A  resolution  cell  is  the  partitioning  of  the  search 
space  in  azimuth  and  elevation.  The  number  of  resolution  cells  that  the  target  spans  is 
predominantly  determined  by  the  system  response  of  the  sensor  system. 

Figure  4.2-1  illustrates  a  typical  IR  system.  The  target  and  background  comprise  a 
The  optical  system  images  the  scene  onto  the  derector  array  which  in  turn  samples  the  intensities 
of  the  image.  The  detector  outputs  are  then  digitized  and  stored  in  a  memory  array  as  the 
digitized  representation  of  the  image.  The  optical  system  is  characterized  by  its  point  spread 
function  (PSF).  This  function  is  the  response  of  the  optics  to  unit  2-D  impulse  functions  or  point 
source  inputs.  Circular  aperture  diffraction  limited  optical  systems  typically  have  an  impulse 
response  with  tire  intensity  in  the  focal  plane  in  a  radial  direction  from  center  given  by: 

E(r)  =  ^[27,  (m)  /  m]2  •  (4.2.1) 


where: 

£0  =  peak  illuminance 

J\  (m)  =  first  order  Bessel  function 

m  a  2je  (NA)  r  /  X 

NA  =  (2  F#)'1  is  the  numerical  aperture  of  the  system 

The  system  aperture  acts  as  a  spatial  filter  with  the  above  characteristic,  bandlimiting  tire 
spatial  frequencies  passed  by  the  optics.  In  essence,  the  pupil  function  that  describes  the  aperture 
is  convolved  with  the  scene  input  to  form  an  image  on  the  detector  array.  The  result  of  the 
spatial  filtering  operation  is  a  reproduction  of  the  scene  with  the  spatial  spectral  content  shaped 
by  the  optics.  In  well  designed  IR  systems,  the  detector  array  makes  use  of  the  lowpass 
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characteristics  of  the  optics  and  samples  the  image  intensities  at  a  spatial  rate  of  at  least  two 
rimes  the  cutoff  frequency  of  the  optics.  This  is  done  to  prevent  aliasing  artifacts.  The  result  is 
an  array  of  Hata  representing  the  image  intensities  as  sampled  by  the  detector  array.  This  is  the 
Hata  that  must  be  analyzed  to  find  targets. 


The  above  IR  system  showed  a  two-dimensional  focal  plane  detector  array  to  simply 
illustrate  the  concept  For  several  reasons  including  detector  technology,  upgrade  potential  and 
cost  many  IR  systems  employ  linear  detector  arrays  and  a  scanning  mechanism  which  scans  the 
image  across  the  array  to  provide  the  second  dimension  to  the  data.  This  is  the  mechanism  used 
in  the  IRAMMP  sensor  described  below.  The  angular  size  of  a  single,  detector  in  image  space 
usually  defines  the  instantaneous  field  of  view  (IFOV)  of  the  sensor  system  and  the  angular 
extent  of  the  image,  scanned  defines  the  total  field  of  view  (TFOV).  The  amount  of  time  that  is 
required  to  cover  the  search  area  with  the  IR  sensor  and  process  the  data  is  known  as  a  scan 
period. 

Because  of  the  nature  of  the  threats,  it  is  sometimes  necessary  to  declare  a  target  within  a 
single  scan  period.  Thus,  this  application  not  only  requires  an  efficient  processing  approach  but 
also  a  robust  detection  performance  in  the  presence  of  clutter  using  only  single  scans  of  data. 
Horizon  IRST  systems  may  not  have  to  contend  with  exceedingly  strong  cloud  and  land  clutter 
or  water  glint.  This  of  course  is  a  function  of  the  extent  of  elevation  coverage  the  sensor 
provides  about  the  horizon.  In  most  cases,  one  may  expect  low  contrast  situations  in  the  vicinity 
of  the  horizon.  In  addition,  IR  clutter  fields  are  in  general  both  spatially  and  temporally 
nonstationary.  This  fact  weighs  heavily  on  any  design  of  a  best  filter  for  detection  processing. 
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The  background  data  that  has  been  used  to  test  filter  algorithms  has  been  extracted  from  the 
Navy's  Infrared  Analysis  Modeling  and  Measurements  Program  (IRAMMP)  database.  The 
IRAMMP  database  contains  an  extensive  collection  of  radiometric  data  including  a  variety  of 
backgrounds.  The  IRAMMP  sensor  used  to  gather  this  data  is  a  dual-band  radiometer  which 
records  IR  backgrounds  in  the  3-5  and  8-12  micron  bands  simultaneously.  For  each  band,  a 
linear  staggered  contiguous  focal  plane  array  consisting  of  120  detector  elements  is  used.  The 
elements  are  staggered  to  allow  alignment  of  the  detector  active  area  edges  in  the  vertical 
direction.  This  allows  for  the  fact  that  if  the  detectors  were  mounted  in  a  vertical  column  with 
their  edges  almost  touching,  there  would  be  a  "dead  space"  between  them  where  no  intensity 
would  be  measured. 

The  detectors  are  electrically  sampled  in  synchronism  with  the  scanning  mechanism  so  that 
successive  samples  of  the  image  taken  from  the  staggered  elements  are  registered  to  the  same 
horizontal  location  in  the  image.  The  instantaneous  field  of  view  (IFOV)  of  each  detector 
element  is  0.23  miiiirariians  vertically  by  0.22  milliradianx  horizontally.  The  impulse  response 
of  the  optical  system  matches  this  closely.  The  optics  were  measured  to  have  85%  of  the  point 
source  throughput  within  1.09  IFOV  for  longwave  and  1.18  IFOV  for  midwave.  The  total  field 
of  view  is  1.6  degrees  vertically  by  5.6  degrees  horizontally.  The  horizontal  scan  is  electrically 
sampled  1480  times  resulting  in  a  spatial  angular  spacing  of  0.0664  milliradians  per  sample. 
This  corresponds  to  a  3.3  times  horizontal  oversampling  of  the  image  by  each  detector  element. 
The  sampling  in  the  vertical  direction  is  not  as  good,  being  only  1  times  (approximately).  The 
result  is  vertical  undersampiing  of  the  image,  with  the  potential  for  spatial  aliasing  artifacts  to 
occur  in  the  vertical  dimension.  This  is  not  a  problem  for  algorithms  that  work  only  on  a  single 
scan  line  basis,  but  the  aliasing  artifacts  could  cause  a  problem  for  2-dimensional  algorithms. 
The  vertical  problem  could  be  alleviated  by  detector  arrays  staggered  to  give  at  least  a  50% 
overlap  of  the  image  as  it  is  scanned  across  the  arrays. 

The  measured  IRAMMP  detector  parameters  are  as  follows:  The  noise  equivalent 
temperature  (NET)  difference  is  0.047  °C  for  the  mid-wave  and  0.032  °C  for  the  long-wave.  The 
noise  equivalent  irradiance  (NED  is  2.6  x  10'4  W/cm2  for  mid-wave  and  2.6  x  10*13  W/cm2  for 
long-wave.  The  dynamic  range  of  the  data  is  84  dB  which  corresponds  to-  15  bits.  Th.. 
uniformity  of  response  across  the  detector  elements  was  not  measured.  It  is  known  that  detector 
point  response  across  the  element  surface  can  vary  quite  a  bit  from  edge  to  edge.  Since  the 
optics  were  well  matched  to  the  detector  size  this  should  not  cause  any  problems.  Inaccurate 
data  results  when  the  optics  create  a  spot  size  much  smaller  than  the  element  size. 
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The  IRAAMP  background  data  can  have  synthesized  or  real  target  data  overlaid  on  it 
digitally  for  testing  various  filter  processes.  Several  different  types  of  detection  filters  have  been 
used  to  date  with  varying  degrees  of  success.  For  the  detection  of  weak  targets  in  homogeneous 
ER  backgrounds,  optimal,  linear,  finite  impulse  response  (FIR)  filters  have  been  developed.  One 
such  filter  is  the  least-mean-squared  filter.  The  nonstationarity  of  IR  clutter  backgrounds 
however,  results  in  degraded  performance  for  linear  filters.  To  overcome  this  deficiency, 
nonlinear  techniques  for  point  target  detection  have  been  proposed.  Within  the  past  several 
years,  a  class  of  nonlinear  signal  processing  algorithms,  collectively  known  as  morphological 
filters,  has  been  applied  to  machine  vision  problems  and,  more  recently,  to  target  detection 
problems.  These  algorithms  respond  to  the  size,  shape  and  orientation  of  imaged  objects  using 
computationally  efficient  logic.  The  result  is  an  improvement  in  real-time  response  and  target 
detection  capability  over  linear  filter  techniques. 

Recent  work  has  shown  that  a  filtering  process  based  upon  a  matrix  formulation  of  the 
Gabor  transform  holds  promise  as  a  mechanism  for  nonlinear  processing  for  the  IRST  problem. 
The  Gabor  transform  can  be  used  as  a  pre-filter  for  the  morphological  processing.  The  resulting 
combination  yields  a  process  that  is  very  much  higher  in  performance  than  linear  algorithms  used 
in  the  IRST  application.  The  next  section  describes  the  formulation  of  the  algorithm  used  and 
details  the  results  of  experiments  with  the  Gabor  transform  algorithm. 

4212  Current  Techniques 

4.2.1.3  Applicability  of  Gabor  Transform 

The  Gabor  transform  is  suggested  for  the  ATR  application  because  it  inherently  resolves  the  data 
subjected  to  it  into  time-  and  frequency  shifted  replicas  of  a  window  function  that  may,  within 
limits,  be  chosen  arbitrarily.  When  the  data  is  presumed  to  include  replicas  of  signals,  a  Gabor 
transform  with  a  high  degree  of  time  resolution  may  be  appropriate  as  an  analysis  mechanism. 
For  unresolved  target  IRST,  the  choice  of  a  window  is  simplified  by  the  fact  that  the  system 
response  to  any  one  target  is  simply  the  point  spread  function  of  the  receive  optics. 

The  optics  of  the  IRST  system  define,  in  general,  the  shape  of  the  response  that  is  expected 
at  the  receiving  end  of  the  system.  Due  to  the  fact  thai  the  system  oversampies  nine  times  for 
each  azimuthal  rotation  of  one  degree,  the  shape  of  the  received  response  to  a  single  point  target 
will  be  well  modeled  by  a  Gaussian  pulse  as  shown  in  Figure  4.2-1. 
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<L2 2  Experimental  Results 

4. 2.2.1  Definition  of  the  Test  Problem 

We  will  now  proceed  to  the  definition  of  the  problem  that  we  used  for  our  Gabor  transform 
experiments.  Fig.  4.2-1  is  a  cloud  scene  from  an  IRST  return  upon  which  we  superimposed  a 
stencil  of  targets  depicted  in  Fig.  4.2-2.  The  resulting  image  is  given  by  Fig.  4.2-3,  in  which  we 
notice  that  the  targets  are  not  discernible  to  the  naked  eye.  This  is  due  to  the  high  dynamic  n 
associated  'vith  this  particular  type  of  scenario.  Next,  we  extracted  one  iine  in  elevation 
contains  about  20  targets  buried  in  the  cloud  cover,  and  used  it  as  a  waveform  (Fig.  4.2-4). 
data  was  then  put  through  a  21-tap  lowpass  filter  and  the  mean  was  subtracted  out,  resulting  in 
the  waveform  of  Fig.  4.2-5,  which  became  the  starting  signal  for  our  1 -dimensional  Gabor 
experiments.  The  question  that  can  now  be  addressed  is  whether  or  not  the  two-dimensional 
time-frequency  representation  that  results  from  the  Gabor  transform  aid  us  in  localizing  the 
targets.  Preliminary  experimental  results  and  detailed  description  of  the  experiments  are 
discussed  in  the  following  sections. 


Fig.  4.2-1.  Goud  cow  section. 


4-27 


m 


Fl«.4^3.  Ckjodcptetmrfctt. 


4.2.Z2  Experiment  Description 

After  importing  the  1024  data  point  signal  described  above  into  the  GSPS  workstation,  we 
proceeded  to  test  different  windows  and  time- frequency  splittings  that  might  yield  results 
deemed  promising  in  either  signal  cleanup  or  target  arrival  detection.  To  this  end  we' began  by 
choosing  the  time-frequency  split  to  be  256  time  points  and  4  frequency  points,  and  used  the 
Longwave  point  response  of  Fig.  4.2-6  as  our  analysis  window.  In  the  two-dimensional  map  of 
Gabor  coefficients.  Fig.  4.2-7,  we  can  clearly  distinguish  isolated  amplitude  peaks  corresponding 
to  some  of  the  target  locations  along  the  off-DC  frequency  lines.  One  will  notice  that  some  of 
the  spikes  have  smaller  amplitudes  than  others  and  some  of  the  targets  are  not  really  well 
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represented  due  to  time  misalignment  of  the  window  with  respect  to  the  signal.  Since  we  have 
chosen  the  number  of  frequency  points  to  be  4,  in  the  critically  sampled  case  it  will  be  sufficient 
to  perform  three  successive  single-point  shifts  of  the  signal  to  guarantee  that  in  at  least  one  of  the 
time  shifts  is  as  well  aligned  as  possible  with  the  chosen  window.  Figs.  4.2-8  through  4.2-10 
represent  the  three  time  shifts,  in  which  we  can  clearly  see  that  different  targets  line  up  better 
than  others  for  a  given  shift  Notice  that  targets  that  initially  lined  up  best  with  the  response 
window  showed  poorly  in  the  shifts.  Each  signal  of  interest  showed  up  most  clearly  in  a  certain 
time  shift. 
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Fig.  4J&-10.  Gabor  eocffldcMs  at  3  shifts. 


In  order  to  address  the  automatic  target  localization  problem,  we  decided  to  perform  some 
simple  nonlinear  processing  to  see  if  we  could  isolate  the  coefficients  that  correspond  to 
target  locations.  This  procedure  was  done  in  the  following  manner.  First  we  decided  to  use 
power  of  MATLAB  to  aid  us  in  this  task,  so  we  saved  the  four  sets  of  coefficients  from  GSPS  in 
MATLAB  format  and  exported  them  into  that  application  to  perform  our  nonlinear  processing. 
Once  we  had  the  coefficients  we  processed  the  coefficients  according  to  the  following  algorithm: 

Let  aO,  al,  a2.  a3  be  the  set  of  Gabor  coefficients  corresponding  to  0,  1,  2,  and  3  time  shifts 
of  the  signal,  respectively.  Then  the  following  pseudo-code  describe  the  process  of  combining 
Gabor  coefficient  values  across  the  shifts: 

for  i=l  to  256 

for  j=  1  to  4 

sum(ja>*  1  a0(j4)  I  + 1  alGa)  I  + 1  a2(j4)  I  + 1  a3(jU)  I 

end 

end 

tnaxvai  =  max  ( sum(j4))  for  all  i.j 
level  =  .05  *  tnaxvai 
for  i  a  l  to  256 
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if  any  ( I  a0(2j)  i ,  i  al(Zi)  I ,  I  a2(2J)  i ,  I  a3{2j)  I  >  lewd ) 
proc_coef(j,i)  *  aO(j4)  for  ail  j 
aid 
aid 

piot(proc_cocf). 


After  the  processing  has  been  done,  import  the  array  proc.coef  of  processed  coefficients  back 
into  the  GSPS  system  and  reconstruct  The  coefficient  map  is  shown  in  Fig.  4.2-11  and  the 
reconstructed  signal  is  shown  in  Fig.  4.2-12.  Notice  that  all  except  one  of  the  targets  were 
isolated,  and  most  of  the  clutter  that  was  originally  in  the  signal  has  been  removed. 
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Clearly  it  will  take  some  more  sophisticated  processing  to  isolate  all  the  targets  and  decrease 
the  false  alarms,  but  as  a  proof  of  concept  this  demonstrates  that  the  Gabor  transform  together 
with  some  nonlinear  processing  techniques  can  be  of  help  in  the  localization  of  targets  in  clutter. 
Of  course  more  work  will  have  to  be  done  to  either  theoretically  or  experimentally  find  the 
optimal  window  as  well  as  tire  optimal  grid  distribution  and  signal  length  in  such  a  way  that  the 
maximal  number  of  targets  is  detected  and  tire  minimal  number  of  false  alarms  is  recorded. 
Other  ways  to  proceed  might  be  to  use  tire  oversampied  Gabor  representation  to  obtain  the  Gabor 
Coefficients,  preprocess  the  signal  further  and  zero  out  all  the  points  that  are  negative,  and  to  use 
the  optimization  methods  described  in  section  4.3  to  find  an  optimal  window  that  is  some  kind  of 
average  of  the  two  window  responses  that  can  be  expected. 

4.2.23  Other  Experiments 

To  ascertain  that  the  above  observations  were  not  artifacts  of  the  uniform  spacing,  or  results 
that  were  somehow  inherent  to  the  cloud  structure,  we  also  ran  the  following  experiments  using 
the  same  analysis  window,  the  same  time- frequency  grid  distribution,  and  the  same  four  time 
shifts  in  order  to  exhaust  the  time  alignment  problem.  We  took  the  same  cloud  sample  as  before 
but  without  any  targets  in  it,  and  used  the  same  analysis  window  to  obtain  Gabor  coefficients  for 
zero  to  three  time  shifts  The  unfiltered  and  filtered  cloud  signals  are  given  in  Figs.  4.2-13  and 
4.2-14,  while  the  time- frequency  maps  are  given  by  Figures  4.2-15-4.2-18.  As  can  be  clearly 
seen,  the  structure  that  is  present  in  the  off-DC  lines  when  the  targets  are  present  is  no  longer 
there  when  the  targets  are  missing,  Le..  the  algorithm  does  not  generate  false  alarms. 
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Fig.  4.2- 14.  Good  filtered. 
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Fig.  40-15.  Gabor  codBriwUa  at  0  shifts. 


Fig.  40-16.  Gabor  cocfTicknls  at  1  shift 
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FI*.  40-17.  Gabor  aMfllcksts  at  2  sklfts. 


The  next  experiment  that  was  performed  was  the  verification  that  the  high  amplitude  Gabor 
Coefficients  were  not  artifacts  of  the  uniform  spacing  of  the  targets. 
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Fif.  4.2-19.  Qoitdi  with  nooHunforiBiy  spaced  taifets. 
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Fig.  4.2-20.  Gabor  coefficients  at  0  shifts. 
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FIf.  4.2-21.  Gabor  cocffiacais  at  1  shift. 
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Flf.  4J2-23.  Gabor  oocfBcMiitsat3  aUfls. 

Even  though  the  window  succeeded  in  pulling  out  the  targets,  the  coefficient  map  is  not  as 
sharp  as  the  map  from  the  uniformly  spaced  coefficients.  We  will  now  see  in  Figs.  4.2-25—4.2- 
28  that  the  use  of  the  Mediumwave  response  window  of  Fig.  4.2-24  will  do  a  much  better  job  of 
pulling  out  the  same  targets. 
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Fig.  42-25.  Gabor  coefficients  at  Osfeifts:  MW 
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Flf.  4J-27.  Gabor  corfflnuti  at  2  shifts:  MW  window. 
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One  can  easily  see  from  the  new  set  of  pictures  that  some  of  the  Gabor  coefficients  that 
correspond  to  target  locations  have  been  magnified,  and  some  of  the  clutter  due  to  spreading  has 
been  eliminated  by  the  use  of  tire  (slightly)  different  window.  Between  the  two  representations, 
all  of  the  targets  that  are  present  in  the  signal  appear  as  large  coefficients  m  the  off- IX!  frequency 
lines.  This  certainly  gives  a  coefficient  map  that  is  very  similar  to  the  one  obtained  from  the 
uniform  target  spacing,  therefore,  we  conclude  that  the  spiky  zones  along  the  off- DC  lines  are  a 
true  manifestation  of  the  presence  of  the  target  and  not  some  artifact  due  to  target  spacing. 

We  also  tried  some  experiments  based  on  the  following  idea.  In  the  impulse  responses 
resulting  from  the  longwave,  the  target  will  always  be  brighter  than  the  clutter  (this  is  not 
necessarily  the  case  with  the  midwave  responses),  and  hence  all  targets  will  have  at  least  a 
portion  which  exhibits  positive  values.  This  implies  that  we  might  be  able  to  take  the  filtered 
signal,  arbitrarily  zero  out  the  points  w!  ich  have  negative  values,  and  process  as  before  with  all 
the  necessary  time  shifts.  The  result  of  this  extra  signal  preprocessing  step  should  be  that  some 
of  the  clutter  will  already  be  removed  before  applying  the  transform,  thereby  cleaning  up  the 
transform  coefficient  map  and  increasing  the  visibility  of  the  signal.  For  all  the  experiments 
described  below  we  used  the  Longwave  response  as  the  window. 

The  result  of  this  zeroing  is  depicted  in  Fig.  4.2-29.  We  then  chose  a  splitting  of  2 
frequency  and  512  time  points  and  obtained  the  Gabor  coefficient  map  for  both  of  the  significant 
time  shifts,  the  result  of  which  is  shown  in  Fig.  4.2-30  (a)  and  (b).  Analysis  using  a  grid  of  4 
frequency  and  256  time  points  was  also  performed.  The  resulting  Gabor  coefficient  map  for  the 
four  relevant  time  positions  is  displayed  in  Fig.  4.2-31  (a)  -  (d).  Comparison  of  the  Gabor 
coefficient  maps  obtained  from  the  preprocessed  and  unpreprocessed  signals  indicates  that  the 
latter  tend  to  spread  the  coefficient  mass  more  across  frequency  and  consequently  cause  the 
target  display  to  become  more  shift  invariant. 
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Fig.  4.2-30.  Gabor  coefficients  for  the  processed  signal  with  M  »  512  and  N  »  2.  (a)  Coefficients  for  tte 
ansMflcd  signal  and  (b)  coefficients  for  the  agw<  after  tint  shifting  once  to  Uw  left. 
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Fig.  4.2-31.  Gabor  coefficients  for  the  procesmd  signal  with  M  ■  256  aad  N  ■  4 .  (a)  -  (d)  correspond  to  the 
coefficients  for  the  signal  after  time  shifting  0  •  3  points  to  the  left  respectively. 
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Fly.  4.2-33.  Oversampled  Gabor  coefficients  for  the  processed  nnshtfted  signal  with  M  *  1024  end  N  =  4. 
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4JL2.4  Results  and  Conclusions 

The  experiments  performed  here  indicate  potential  for  Gabor  analysis  as  a  component  of 
unresolved  target  detection  for  the  IRST  problem.  Synthetic  targets  representing  replicas  of  the 
receive  optics  point  spread  function  were  distributed  both  uniformly  and  nonuniformly  across  an 
observed  cloud  background,  and  were  found  to  be  visible  in  processed  Gabor  expansions  that 
favored  time  resolution  as  opposed  to  frequency  resolution.  There  is  a  robustness  in  the  results 
that  indicates  stability  of  the  findings  with  respect  to  the  details  of  the  approach. 
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<L2J3  Future  Work 

Extensions  of  the  current  work  would  emphasize  several  things,  among  them  further  theoretical 
studies  to  gain  an  explanation  of  the  moderate  successes  observed  in  the  reported  experiments. 
Because  of  the  time  resolution  inherent  in  the  oversampled  approach,  further  experiments  with  it 
are  warranted  also.  Application  of  the  techniques  to  additional  data  sets  could  be  expected  to 
^nhanft.  our  understanding  of  the  performance  as  well  as  uncover  areas  where  the  techniques 
shown  here  are  either  inadequate  or  in  need  of  support  from  other  methodologies.  In  particular, 
we  are  anxious  to  investigate  algorithms  which  combine  the  power  of  Gabor  methodologies  in 
producing  fa™  that  exhibits  signal  detection  in  a  few  large  coefficients  and  morphological 
filtering  to  aggregate  the  observed  points  and  replace  some  “human  operator”  functions  with 
machine  capability. 
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4J  MINIMUM  DIMENSION  GABOR 

Under  contract  to  DARPA,  AA EC  earned  out  Phase  I  and  Q  SBIR  studies  on  the  topic  of 
“Minimum  Dimension  Gabor  Representations.”  This  work  dnves  at  obtaining  Gabor  expansion* 
of  a  set  of  signals  in  which  the  majority  of  the  information  about  the  signals  is  concentrated  in  a 
few  large  expansion  coefficients.  Methods  of  nonlinear  optimization  theory  are  employed  in  this 
quest  The  work  carried  out  under  these  contracts  produced  very  positive  results  leading  to  two 
conference  publications  [SWEE,  ORR],  but  did  not  complete  the  needed  research.  Some  of  the 
research  extensions  were  carried  out  under  this  contract  and  their  results  are  presented  below. 

In  the  following  paragraphs  we  summarize  the  problem  statement  the  theoretical  results 
produced  under  the  Phase  I  SBIR  and  the  numerical  algorithms  and  results  from  Phase  II. 
Following  this  the  new  work  performed  under  this  contract  is  described. 

4.3.1  Summary  of  SBIR  Contract  Effort 

43.1.1  Problem  Description 

Series  expansions  of  signals  in  which  significant  features  of  the  signal  are  captured  in  a  few  large 
coefficients  are  desirable.  This  work  shows  that  given  a  collection  of  signals,  it  is  possible  to 
find  Gabor  representations  for  these  that  are  maximally  concentrated  in  time-frequency  space. 
The  problem  addressed  is:  given  a  signal  set,  find  the  window  function  of  the  Gabor  expansion 
that  minimizes  an  “average  dimension”  of  the  signal  representations  relative  to  that  window. 
The  dimension  measure  employed  is  entropy  based  and  related  to  the  quantum-mechanical 
technique  where  one  interprets  expansion  coefficients  as  probabilities. 

An  iterative  algorithm  based  on  partial  derivatives  of  the  signal  set  dimension  with  respect 
to  the  expansion  function  was  used  to  evaluate  effectiveness  of  several  nonlinear  optimization 
algorithms  in  finding  an  optimum  window. 

4 .3.1.2  Phase  I  Accomplishments 

The  key  results  of  Phase  I  are  presented  here,  and  the  reader  is  referred  to  the  final  report 
and  a  conference  publication  for  details  [AAEC9l(2),  ORR].  Key  findings  include: 

*  Proof  that  the  notion  of  "dimension"  is,  for  a  single  signal,  completely  arbitrary,  depending 
wholly  upon  the  choice  of  representation  and  not  at  ail  on  the  signal; 

*  A  convincing  demonstration  that  where  a  set  of  signals  is  concerned,  the  structure  of  the  set 
bears  an  inherent  relation  to  a  dimension  that  can  be  assigned — this  dimension  remains  a 
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function  of  the  representation,  but  can  be  over  quite  general  subsets  of  all  basis 

functions; 

•  Recognition  that  the  constrained  structure  of  tune-frequency  basis  sets  such  as  the  Gabor, 
affine  wavelet,  etc.,  make  these  ideal  classes  of  bases  over  which  to  carry  out  the 
minimization  of  set  dimension; 

•  A  determination  of  a  reasonable  set  of  requirements  for  the  definitions  of  both  signal  and 
signal  set  dimensions; 

•  A  derivation  from  these  requirements  of  a  signal  dimension  formula  that  is  unique  to  within 
a  single  parameter  and  a  set  dimension  formula  that  is  then  unique  to  within  the  assignment 
of  a  weighting  function — the  formula  combines  the  quantum  mechanical  relationship  of 
expansion  coefficients  to  probabilities  and  the  information  theoretic  notion  of  the  entropy 
of  a  probability  distribution: 

•  Successful  implementation  of  the  signal  dimension  measure  within  AAEC's  Gabor 
processing  software  testbed  (GSPS); 

•  A  useful,  though  incomplete,  classification  of  signal  set  types  in  a  manner  that  motivates 
physically  meaningful  choices  of  weighting  functions  for  dimension  assignment; 

•  An  extensive  compilation  of  examples  of  basis  set  constraints  that  would  allow 
optimization  over  more  restricted  subsets  of  Gabor  or  other  representation  families: 

•  Development — via  use  of  Poisson  summation  and  the  ambiguity  function— of  new 
characterizations  of  Weyl-Heisenberg  frames  (the  immediate  generalization  of  Gabor 
representations)  that  should  make  it  easier  to  perform  dimension  optimization  under 
constraints  within  these  structures. 

These  achievements  set  the  starting  point  for  Phase  II.  Having  successfully  formulated  the 
dimension  concept  and  the  associated  minimization  problem,  the  immediate  need  was. to  discover 
the  extent  to  which  analytic  machinery  can  be  brought  to  bear  on  solving  the  difficult  nonlinear 
optimization  problem  we  have  created.  For  example,  in  basis  systems  that  are  defined  by  a 
single  window  function,  which  include  most  of  the  cases  of  interest— Gabor,  Weyl-Heisenberg, 
affine  wavelets,  etc. — the  solution  consists  of  determining  an  optimum  window.  If  there  are 
sufficiently  restrictive  constraints  placed  upon  the  class  of  eligible  windows,  it  should  be 
possible  to  obtain  solutions  that  not  only  minimize  the  cost  function,  but  have  other  desirable 
behaviors  as  well.  If  we  can  obtain  analytic  answers  for  some  simple  but  nontrivial  cases,  the 
utility  of  the  techniques  will  be  enhanced.  To  conclude  this  aspect  it  is  important  to  achieve  a 
software  capability  that  implements  the  optimization  algorithms  that  are  found. 
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Phase  Q  concentrated  mainly  on  application  of  a  concept  put  forth  under  the  prior  phase 
[AAEC92(2)j.  New  theoretical  developments  woe  limited  to  only  those  needed  to  cany  out  the 
optimizations  or  to  interpret  results.  The  main  theoretical  idea  is  of  course  the  development  of 
the  discrete  domain  cost  function  for  optimization  and  the  choice  of  optimization  variables. 
Carrying  out  the  optimizations  over  the  biorthogonal  function  has  proven  quite  successful,  and 
has  ni«Tnnim»ri  an  intermediate  computational  burden,  one  that  could  have  been  quite  expensive, 
were  the  window  itself  directly  optimized.  We  see  nothing  in  the  results  to  prompt  any  revision 
our  underlying  methodology. 

The  final  theory  topic  is  the  evaluation  of  computational  complexity  of  the  dimension 
functions.  We  found  that  a  dimension  or  gradient  evaluation  for  one  signal  is  an  OfP2)  process 
(recall  that  P  is  the  number  of  points  in  a  signal),  whereas  a  Hessian  computation  is  0(P3). 
These  are  reasonably  high  orders — think  of  the  effort  people  will  expend  to  replace  a  traditional 
order  Of/*2)  Fourier  transform  with  an  FFT.  No  pressing  need  to  rely  on  great  numbers  of 
Hessians  has  been  uncovered,  so  we  can  probably  say  with  reasonable  confidence  that  the  per- 
signai  computation  caps  at  0(Pa).  No  analysis  of  the  adaptive  algorithms  was  attempted  since 
there  are  so  many  variants,  and  the  run  time  is  data-dependent. 

Initial  experiments  using  a  local  gradient  method  were  very  successful  in  demonstrating  that 
iterative  techniques  can  converge  to  a  good  solution.  A  number  of  other  things  were  observed 
from  these  as  well.  Convergence  to  a  suboptimum  point,  i.e.,  a  local  minimum,  was  seen 
occasionally,  and  in  some  circumstances  the  answer  was  a  very  plausible  candidate.  In  some 
cases  the  converged  biorthogonal  caused  the  Gabor  coefficients  of  the  signal  to  be  dominated  by 
a  single  large  coefficient  located  somewhere  other  than  the  origin;  the  resulting  window  was 
either  incalculable  or  did  not  resemble  the  signal.  This  circumstance,  should  be  regarded  as  a 
successful  intermediate  result  from  which  the  search  for  a  more  properly  behaved  solution  can  be 
sought.  Only  partial  success  in  carrying  on  from  this  point  was  achieved,  leaving  this 
circumstance  as  one  of  the  candidates  for  further  work. 

The  more  appropriate  experiments  using  the  optimization  techniques  of  Matlab  and  NAG 
yield  further  confirmation  that  the  process  would  work,  but  some  more  surprising  results  were 
found.  Most  startling  among  these  is  the  observed  convergence  to  highly  unstable  answers. 
There  exist  cases  for  which  a  solution  of  dimension  1.0  is  found  using  a  window  having  a  very 
large  dynamic  range  of  values,  on  the  order  of  1015  in  some  cases.  In  general  these  windows 
tend  to  pile  up  at  the  end  of  the  interval,  and  yield  Gabor  coefficient  distributions  in  which  the 
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large  coefficient  is  offset  to  the  last  time  point  but  not  offset  in  frequency.  The  instability  of 
these  is  seen  when  the  attempt  to  reconstruct  the  signal  in  GSPS  after  clipping  off  the  low 
coefficients  fails  violently,  apparently  due  to  numerical  errors.  The  NAG  routines  usually  reach 
a  happy  termination  in  these  cases,  confirming  that  the  cost  function  bong  used  is  insensitive  to 
the  behavior. 

In  other  cases,  the  optimizations  terminate  at  low  (well  under  2.0)  dimensions  and  find 
windows  that  resemble  the  signal  except  for  some  usually  spiky  irregularities.  These  runs  seem 
to  be  finding  local,  not  global,  minima,  and  are  cause  for  some  concern,  because  the  excessive 
behaviors  of  some  of  the  other  undesired  termination  types  are  not  present  here.  On  the  other 
hand,  the  results  are  not  unusable.  Convergence  to  slightly  suboptimum  answers  may  in  many 
cases  be  almost  as  good  as  finding  the  best  answer.  Although  we  would  be  well  advised  to 
understand  the  phenomema  at  work  here,  it  is  not  clear  to  what  extent  these  results  represent 
problems  to  solve  in  the  future. 

The  issues  raised  by  the  single  signal  cases  remain  when  multiple  signals  are  present.  There 
are  cases  in  which  the  window  ‘blows  up’  as  described  earlier,  but  still  yields  a  low  dimension 
representation  for  each  signal. 

» 

In  other  respects  the  results  are  encouraging.  In  this  light  we  cite  the  experiments  in  which 
the  optimum  window  was  sought  for  a  signal  set  containing  time  translates  of  a  decaying 
exponential.  The  dimension  of  such  sets  was  computed  for  families  of  exponential  windows, 
and  the  intuition  about  those  results  is  that  the  exponential  window  is  not  quite  optimum.  It  was 
not  until  the  multi-signal  optimization  code  was  working  that  we  could  look  for  the  optimum  in 
such  a  case.  When  we  did,  we  found  a  window  with  strong  exponential  characteristics  that  stably 
represented  all  eight  signals.  The  dimension  of  the  set  was  similar  that  found  with  the  best 
strictly  exponential  window. 

New  issues  encountered  during  the  multi-signal  runs  relate  to  computation  time.  The  effort 
•  to  compute  dimension  or  its  derivatives  scales  linearly  with  the  size  of  the  signal  set.  and  there  is 
no  apparent  way  to  decrease  this.  Overall  run  time  is  somewhat  unpredictable  despite  our 
understanding  of  the  computational  burden  of  dimension-related  evaluations,  because  the  number 
of  evaluation  rails  is  both  data-  and  method-dependent.  These  observations  tend  to  confirm  the 
suspicion  that  optimum  window  analysis  is  an  off-line  activity.  This  is  not  at  all  a  poor  finding, 
given  that  many  of  the  envisioned  applications  of  the  techniques — e.g.,  signal  detection  and 
classification — might  fall  into  such  categories. 
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Results  of  the  Mailab  minimization  nuts  indicate  that  the  choice  of  algorithm  is  somewhat 
H«ta  dependent.  From  the  six  algorithms  tested,  two  minimized  the  dimension  in  all  the 
experiments,  and  the  other  four  algorithms  failed  for  certain  experiments.  Therefore  the  first  two 
algorithms  can  be  considered  more  reliable,  but  not  necessarily  the  optimum  in  all  cases. 

It  was  also  observed  that  the  number  of  iterations  was  reduced  by  using  the  gradient  if  it  was 
available.  When  the  gradient  was  supplied,  the  Cubic  Interpolation  line  search  method  was  seen 
to  perform  better  than  Mixed  Polynomial  Interpolation;  fewer  function  evaluations  were 
performed,  but  more  gradients  were  evaluated. 

The  Simplex  Search  minimization  method  performed  the  worst,  resulting  in  a  successful 
minimization  in  only  two  out  of  six  experiments.  In  the  two  where  it  did  finish,  the  number  of 
iterations  was  the  largest  of  all  the  methods. 

The  NAG  optimization  routines  have  been  shown  to  be  capable  of  determining  Gabor 
representations  with  minimum  dimension.  This  capability  has  been  demonstrated  both  for 
simple  functions  and  more  complex  real  world  signals,  and  for  single  and  multiple  signals. 

The  initial  estimate  of  the  biorthogonal  is  only  marginally  important  in  controlling  if.  and 
how  efficiently,  the  optimizing  routine  finds  a  solution. 

In  general,  the  quasi-Newton  optimization  method  employing  a  user  supplied  first  derivative 
of  the  cost  function  (dimension)  yielded  the  best  results.  It  required  fewer  iterations  to  converge 
than  the  quasi-Newton  method  with  the  finite  difference  approximation  to  the  first  derivative, 
and  ran  considerably  faster  than  the  Newton  method  which  employed  a  user  supplied  second 
derivative.  Each  method  seemed  equally  capable  of  finding  solutions. 

Linear  constraints  on  the  cost  function  during  optimization  were  helpful  when  a  priori 
information  about  tire  biorthogonal  was  available.  For  the  more  complex  transient  signals,  where 
the  optimum  biorthogonal  is  not  known,  linear  constraints  were  not  helpful. 

The  best  number  of  frequency  and  time  points  (M  and  AO  for  a  given  number  of  samples 
depends  somewhat  on  the  signal.  For  signals  with  detailed  features  near  the  beginning,  better 
solutions  were  generally  obtained  when  M>  N.  For  signals  without  meaningful  content  near  the 
beginning,  Af  <  N  generally  produced  better  results. 

Solutions  were  found  for  the  multiple  signal  experiments  where  the  average  dimension  was 
less  than  if  one  of  the  signals  were  employed  as  the  window  function. 
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Use  of  the  free  parameter  in  the  dimension  formula  (exponent  r )  was  not  fully  exploited  in 
the  tests,  but  we  were  able  to  determine  some  consequences  of  varying  it.  Larger  values  of  r 
were  found  to  speed  convergence  in  cases  that  were  already  convergent  at  smaller  r.  Large  r 
tends  to  amphasige  the  greatest  Gabor  coefficient  and  push  the  distribution  towards  the  desired 
consequently,  the  optimizer  sometimes  behaves  as  though  its  job  has  become  easier,  and 
presents  results  that  appear  to  have  stopped  short  of  reaching  the  ‘eyeball  optimum.  *  The  overall 
conclusion  is  that  low  r  is  the  most  sensitive  case,  and  may  be  best  for  fine  tuning  a  result,  while 
a  larger  value  of  r  might  be  used  to  get  to  the  vicinity  of  an  acceptable  solution.  Algorithms  that 
automatically  adjust  r  en  route  according  to  criteria  related  to  these  observations  can  be  imagined 
for  the  future. 

4.3.13  Conclusions 

We  have  shown  that  it  is  possible  to  begin  with  a  collection  of  signals  and  find  a  Gabor 
representation  of  these  that  is  maximally  concentrated  in  the  sense  of  the  dimension  function 
H*tin»ri  in  Phase  I.  Some  of  the  computed  examples  show  significant  differences  between  a 
window  function  found  using  optimization  techniques  and  a  more  naively  chosen  one,  some  do 
not.  As  expected,  the  techniques  developed  here  are  not  automatically  candidates  for  every 
application.  Tn««*d.  they  provide  a  body  of  technique  that  may  enhance  the  ability  to  carry  out  a 
few  procedures.  Further  evaluation  must  await  their  application  to  some  tasks. 

4.3.2  Extensions  Under  Current  Contract 

43.2.1  Optimization  Tools 

The  extensions  under  this  contract  were  a  continuation  of  the  nonlinear  optimization  experiments 
using  a  more  complex  and  powerful  NAG  routine,  E04UCF  than  the  previously  used  algorithm 
E04KAF.  The  important  properties  of  the  E04KAF  routine  are  displayed  below: 

•  Quasi-Newton  algorithm 

•  Uses  analytic  gradient  (first  derivative) 

•  Builds  up  surface  curvature  information  (Hessian,  or  second  derivative) 

•  Incorporates  bounds  on  independent  variables 

The  NAG  routine.  E04UCF,  is  more  powerful  because  it  allows  the  user  to  control  more 
optimization  parameters.  It  performs  nonlinear  optimization,  including  bounds  on  the  variables, 
linear  constraints,  bounds  on  the  linear  constraints,  nonlinear  constraints,  and  bounds  on  the 
nonlinear  constraints.  It  uses  a  sequenuai  quadratic  programming  algorithm  in  which  the  search 
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direction  is  the  solution  of  a  quadratic  programming  problem  (NAG).  Its  properties  are  listed 
below: 

•  Sequential  quadratic  programming  algorithm 

•  Minimizes  smooth  nonlinear  function  subject  to  constraints 

•  Takes  analytically  specified  first  pardals 

•  Approximates  unspecified  first  pardals  by  finite  differences 

•  Incorporates  linear  and  nonlinear  constraints 

•  Incorporates  bounds  on  independent  variables,  linear  constraints,  nonlinear  o 

A3JJ2  Experiments 

Numerous  experiments  were  performed  with  the  E04UCF  optimization  routine  to  determine 
average  windows  for  several  test  cases. 

The  NAG  routine  E04UCF  was  used  with  several  combinations  of  bounds,  linear  and 
nonlinear  constraints.  Experiments  were  performed  with  two  signal  sets:  a  single  signal,  and  two 
signals.  The  one  signal  case  used  a  real  signal  chosen  from  the  first  data  set  of  the  fault 
identification  section,  with  a  Gabor  lathee  of  Af  =  4,  N*  16.  The  two  signal  case  was  composed 
of  a  rectangular  pulse  and  a  decaying  exponential,  with  M  »  4,  N  »  4. 

The  bounds  and  constraints  were  chosen  from  the  available  set  of: 

•  Bounds  on  the  variables 

•  Linear  constraints 

•  Bounds  on  linear  constraints 

•  Nonlinear  constraints 

•  Bounds  on  nonlinear  constraints 

The  choice  of  these  bounds  and  constraints  permitted  the  control  of  the  minimization  process  to 
avoid  convergence  to  unstable  results. 

Three  types  of  constraints  were  used  and  they  are  described  below: 
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•  Type  1:  Nonlinear.  Restricts  exponential  growth  of  bionhogonal  signal  envelope 

Desired  constraint  is  of  the  form  shown  in  Fig.  4.3-1.  but  hard  to  use  because  of  lack  of 
smoothness 


s  amax{|£5j,|fc6|)|k7j,|£>8|} 
max{|i>5|,|i>6|,|i>7|,|i>8|}  S  amax{|d9|,|610j)|fcu|,|Z>12|} 
max{|&9|,|&10j,jiJ11|,|£12|}  >  amax{|i13|,|*14|,j*15|,|i»16j} 


Fig.  43-1  Nonlinear  constraint  that  controls  exponential  growth  of  the  biorthogonal  on  the  basis  of  the 
largest  magnitude  value  within  each  Gabor  time  slice. 
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This  type  of  constraint  was  replaced  by  pointwise  comparisons  as  shown  in  Fig.  43-2; 


M> 

N> 

kg  |  > 


«P5|. 

a\b9\, 


1*2 1  >  a|*s|.  I*sl  >  «l*r|.  h\  >  afel 

*6 1  >  «l*lol-  l*7l  >  al*lll>  1*8 1  ^{*12! 

*lo|  >  a|*i4|,  I  i’ll  I  >  ^1*15 1>  1*12!  >  <*1*16 


Fig.  4-3-2  Nonlinear  constraint  am  control  exponential  growth  of  the  biorthogonal  on  the  bee u  of  the 
nmnitedM  point-by-point  within  cacti  Gabor  time  slice. 
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•  Type  2:  Nonlinear.  Controls  dynamic  range  of  biorthogonal  within  first  Gabor  time 
interval.  The  need  for  this  can  be  seen  by  examining  the  inversion  that  the  window 
function  from  the  biorthogonal  in  the  matrix  Gabor  method.  Not  usable  due  to  non¬ 
smoothness.  Not  needed  thus  far. 


Fig.  43-3.  Nonlinear  constraint  that  controls  range  within  tfae  first  Gabor  time  interraL 
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Type  3:  Linear  Restrict  positive  exponential  growth  of  biorthogonal 


>  ai\ 3,  &2  >  od> 14,  b$  >  a^j5,  b4  >  cxbl6 


Fig.  4-3-4.  Linear  conttraiiit  tot  control  wpannitM  growth  of  the  btorthofoml  on  the  b«i«  of  the  lirpit 
magnitude  nine  within  the  first  and  last  Gabor  time  dice. 

The  objective,  or  cost  function,  computed  the  minimum  dimension  of  the  signal  set  for  each  of 
the  experiments  shown  in  the  following  figures.  For  the  single  signal  figures,  the  initial 
biorthogonal  was  a  rectangular  pulse. 

The  single  signal  cases  were  made  with  a  transient  signal  from  the  first  data  set.  Fig.  4.3-5 
shows  the  constraints  used,  the  signal,  minimized  biorthogonal,  the  window,  the  Gabor 
coefficients,  and  the  dimension  was  found  to  be  1.06.  This  appears  to  be  an  excellent 
minimization  result,  because  the  dimension  is  very  close  to  1.  However,  the  resulting  window  is 
a  single  negative  spike  with  amplitude  of  1.7e+21,  not  a  useful  window.  This  is  the  result  of 
convergence  to  an  unstable  minimum.  This  indicates  that  the  minimization  process  needs  to  be 
constrained  to  avoid  convergence  to  this  kind  of  unstable  minimum 

Fig.  4.3-6  shows  the  same  case  as  in  the  previous  figure',  except  that  the  NAG  routine  E04UCF 
was  used,  and  nonlinear  bounds  and  constraints  were  applied.  The  resultant  window  was 
reasonable,  and  the  dimension  was  still  low,  at  1.75. 
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Single  signal  -  Method  E04KAF 


Bounds  on  variables 


Linear  constraints  _  _  N 


Bounds  on  linear  constraints _  _  N 


Nonlinear  constraints  N 


Bounds  on  nonlinear  constraints  N 


Bounds  on  nonlinear  constraints  N 


D  ~  1.06 


Fig.  43-5.  Transient  signal  -  bounds,  no  constraints 
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Single  signal  •  Method  E04UCF 


D  =  1.75 


4-62 


Rna I  TR:APPUCA*nON  OF  THE  GABOR  REPRESENTATION  TO  MHJTARY  PROBLEMS 


Several  expenments  were  performed  with  a  two  signal  set:  a  rectangular  pulse,  and  a 
delayed  exponential.  For  all  of  the  two-signal  minimization  experiments,  a  random  initial 
biorthogonai  was  chosen.  This  initial  biorthogonal  and  its  corresponding  window  are  shown  in 
Fig.  4.3-7. 


Before  attempting  to  find  an  optimum  window,  each  of  the  two  signals  were  used  as  a 
window  in  the  computation  of  Gabor  coefficients,  as  shown  in  Fig.  4.3-8  and  4.3-9.  In  Fig.  4.3- 
8,  the  pulse  is  used  as  a  window,  which  happens  to  be  the  same  as  the  biorthogonal.  The  lower 
coefficient  set  reproduces  the  pulse,  and  the  upper  set  reproduces  the  exponential.  The  average 
of  the  two  dimensions  is  1.98.  Similarly,  in  Fig.  4.3-9,  the  exponential  used  as  a  window  results 
in  a  average  dimension  of  1.303.  What  we  are  looking  for  is  an  optimum  window  that  produces 
a  minimum  average  dimension.  Figure  4.3-10  shows  the  optimization  result  using  the  NAG 
routine  E04KAF,  which  uses  bounds  on  the  variables  only.  The  resulting  average  dimension  is 
1.496,  very  close  to  the  previous  case. 

The  next  experiment  (Fig.  4.3-1 1)  used  the  E04UCF  method  with  bounds  on  the  variables 
and  linear  constraints  with  bounds.  The  four  linear  constraints  used  were  the  Type  3  constraints 
mentioned  earlier.  The  result  was  D  =  2.00,  showing  that  this  type  of  constraint  did  not  succeed 
in  decreasing  the  average  dimension. 

Figure  4.3-12  shows  the  result  of  using  bounds  on  the  variables  and  twelve  Type  1  nonlinear 
constraints  with  bounds.  The  average  dimension  is  lower,  at  1.887. 

The  value  of  at,  the  coefficient  in  the  Type  1  constraint  equation  was  then  changed  from  0.5 
to  1.01.  The  previous  case  was  run  again,  resulting  in  Fig.  4.3-13.  The  minimized  biorthogonal 
was  constrained  as  expected,  but  the  average  dimension  was  higher.  The  same  case  was  run 
again  with  alpha  =  1.02,  with  the  result  shown  in  Fig.  4.3-14.  The  average  dimension  is  now 
lower,  but  the  biorthogonal  has  changed  to  nearly  a  rectangular  pulse.  Apparently,  the 
optimization  algorithm  is  very  sensitive  to  changes  in  the  constraint  function. 
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Imtiai  biorthogonal  chosen  randomly  (upper  panel) 
Corresponding  window  (lower  panel) 


F!g.4J»7.  Initial  Morthofo—l  and  conopondlog  window. 


4-64 


nnat  TRlAPPUCATION  OP  THE  GABOR  REPRESENTATION  TO  IfSUTARY 


Two  signals  -  No  optimization 
Window  equals  rectangular  pulse 
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Two  signals  -  No  optimization 
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Two  signals  -  Method  E04KAF 


Bounds  on  variables 

■ 

Linear  constraints 

N 

Bounds  on  linear  constraints 

N 

Nonlinear  constraints 

N 

Bounds  on  nonlinear  constraints 

N 

Two  signals  -  Method  E04UCF 


Bounds  oo  variables 


Linear 


Bounds  on  linear  constraints 


Nonlinear  constraints 


Bounds  on  nonlinear  constraints 
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Two  signals  -  Method  E04UCF  (a*J) 
- 12  constraints 


Bounds  on  variables 


Linear  constraints 


Bounds  on  linear  constraints 


Nonlinear  constraints 


Bounds  on  nonlinear  constraints 
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Fig.  4.3-12.  Two  signal  case  •  bounds,  nonlinear  constraints. 
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Two  signals  -  Method  E04UCF  (a  *  1.01) 


Dw  =  2.35 


Fif.  43.13.  Two  signal  cate  •  bonds,  nonlinear  constraints. 
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Two  signals  -  Method  E04UCF  (a  =  1.02) 


Bounds  on  variables 


Linear 


Bounds  on  linear  constraints 


Nonlinear  constraints 


Bounds  on  nonlinear  constraints 
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Dav=  1.97 

Fit.  43-14.  Two  «ana!  case  •  bonds.  non 
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With  these  additional  optimization  runs,  we  have  managed  to  minimize  the  dimension  where  the 
biorthogonal  function  previously  did  not  converge  to  a  reasonable  state.  The  addition  of 
nonlinear  constraints  and  bounds,  along  with  the  sequential  quadratic  programming  minimization 
algorithm,  has  given  us  much  more  control  over  the  optimization  process.  It  was  observed  that 
the  minimization  can  be  very  sensitive  to  changes  in  the  constraints.  This  was  seen  in  the 
significant  change  in  the  resulting  window  function  resulting  from  a  small  change  in  the  value  of 
alpha  in  the  nonlinear  constraint  equation. 

4.3.3  Future  Work 

Tire  work  presented  above  makes  a  positive  extension  of  the  results  originally  obtained  under  the 
SBIR  and  provides  further  justification  of  the  ideas  that  first  lead  to  this  line  of  endeavor.  The 
task  remains  unfinished,  however.  Research  to  date  has  not  completed  the  task  of  determining 
the  best  constraints  to  use  in  optimizing  windows,  and  the  bugs  in  using  the  optimization 
software  have  not  yet  been  completely  worked  out  Work  has  been  slowed  somewhat  by  this 
situation.  Using  simulated  annealing  for  optimization  has  also  been  suggested,  and  AAEC  is 
currently  looking  into  available  software  packages  for  that  purpose. 

In  summary,  the  effort  so  far  has  proved  in  principle  most  of  the  supporting  concepts,  but 
has  been  insufficient  to  transition  the  work  into  the  applications  arena  as  yet.  AAEC  sees 
particular  promise  for  this  technology  in  certain  applications  areas,  and  is  planning  to  propose 
effort  in  those  areas.  A  key  area  is  automatic  target  recognition  (ATR).  Initial  work  in  applying 
the  Gabor  transform  to  problems  within  that  discipline  is  reported  elsewhere  in  this  document. 
Machine-aided  recognition  problems  have  the  feature  that  searching  for  objects  can  be  enhanced 
in  circumstances  where  shape  characteristics  of  the  objects  are  partially  known  in  advance,  either 
through  a  priori  knowledge  or  data-aided  algorithms.  For  example,  in  signal  analysis,  the  Gabor 
transform  is  particularly  adept  at  finding  features  having  a  common  envelope. 

To  maximally  exploit  such  a  circumstance,  extraction -of  the  analysis  window  from  the  data 
looms  important.  Given  a  large  body  of  data  such  as  that  often  encountered  in  an  ATR  problem, 
use  of  the  data  to  drive  the  analysis  functions  seems  wise  as  a  measure  to  cut  the  amount  of  blind 
search,  especially  in  view  of  findings  that  allegedly  more  “robust”  tools  such  as  the  Wigner 
distribution  can  create  artifacts  through  nonlinear  processing  if  not  used  carefully.  The  role  for 
optimum  Gabor  windowing  in  this  scheme  is  clear,  and  as  a  result  it  appears  that  the  best  way  in 
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which  to  continue  the  line  of  wok  discussed  above  is  to  do  it  within  the  context  of  an  application 
area  such  as  ATR. 

The  research  is  at  the  point  where  it  could  profit  from  the  interaction  with  real  data  as  an  aid 
in  algorithm  development/refinement.  AAEC  anticipates  proposing  a  body  of  work  of  this 
nature  as  a  logical  follow-on  the  work  performed  in  both  this  contract  and  the  cited  SBIR’s. 
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STATEMENT  OF  WORK 

AAEC  will  perform  tbe  following  work  items,  segregated  by  year. 

Year  1: 

•  Proposed  Theoretical  Developments:  Study  methods  to  improve  the  stability  and 
accuracy  of  the  Gabor  coefficient  computations;  perform  error  analyses  to  understand  tbe 
accuracy  of  approximation  achieved  by  truncated,  finite  Gabor  expansions. 

•  Software  Development:  Develop  software  code  on  a  GFE  Aspen  computer  to  execute 
the  new  algorithms  for  Gabor  analysis.  This  includes  full  debugging  and  achieving  operational 
status. 

•  Comparative  Analysis  of  Competing  Methods:  Investigate  tbe  utility  of  competing 
methods,  including  other  double  series  representations.  Wigner  distributions,  non-abeiian 
harmonic  analysis,  etc.,  for  the  problems  under  attack. 

•  Potential  Applications:  Begin  to  investigate  potential  application  of  Gabor 
representations  to  military  problems,  and  use  the  tentative  results  to  guide  the  needed 
theoretical  developments.  Test  problems  to  be  chosen  from  those  in  the  proposal 

Year  2: 

•  Potential  Applications:  Continue  to  investigate  potential  application  of  Gabor 
representations  to  military  problems,  and  use  the  tentative  results  to  guide  the  needed 
theoretical  developments.  Test  problems  to  be  chosen  from  new  areas  discovered  during 
first-year  investigations  and  discussions  with  interested  parties  at  DARPA  and  elsewhere. 

•  Software  Upgrades:  Maintain  existing  software  and  add  new  capabilities  as  needed  to 
improve  ability  to  model  potential  applications. 


Year  3: 

•  Continuation  of  Year  2  work  items. 

Full  Statement  of  Work  disclosed  to  DARPA  and  available  upon  request. 
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COMPUTATIONAL  ACCURACY  AND  STABILITY  ISSUES 
FOR  THE  FINITE.  DISCRETE  GABOR  TRANSFORM 

Rogtlio  Baton  and  Richard  S.  Orr 

Atlantic  Aerospace  Electronics  Corporaooa.  6404  Ivy  Lane.  Suite  300.  GreenbdL  MD  20770  USA 


ABSTRACT 

A  Gabor  expansion  may  employ  highly  nonorthogonai 
basts  functions  and  consequently  indent  accuracy  and 
stability  problems  due  to  near-singulamy  when  computed 
digitally.  Two  methods  to  discretize  the  Gabor  transform 
are  studied  from  tne  viewpoint  of  controlling  numerical 
properties:  (i)  a  Zak  transform-based  method  and  (itt  a 
matrix  method.  Theoretical  issues  relating  to  the  singular 
behavior  of  each  arc  cued,  and  stabilization  techniques  are 
proposed.  We  demonstrate  the  validity  of  each  technique 
with  results  of  numerical  experiments,  concluding  that 
stability  and  accuracy  can  usually  be  achieved  in  a 
digitally  implemented  Gabor  transform  by  proper  choice 
of  algorithm  and  stabilization  mechanism. 

L  INTRODUCTION 

A  Gabor  expansion  of  an  LHK)  function  is  a  doubly 
infinite  senes  of  the  form 

/(')»  1  Z**u,»nu,W 


in  which  the  are  the  Gabor  coefficients,  and  the 

lw-.l  are  time  and  frequency  translations  of  a  window 
function,  wit),  that  are  the  basis  of  the  expansion.  The 
properties  of  the  expansion  are  highly  dependent  upon  the 
behavior  of  the  window  with  respect  to  the  ttme- 
frequency  lattice  » [(mIT.nT))  that  supports 

the  basts.  The  parameter  T  >  0  mitigates  the  tune  and 
frequency  resolution  of  the  transform.  Many  of  these 
properties  are  readily  expressed  in  terms  of  the  Zak 
transform  of  w. 


Zwiv.r)*  2wf*7*+T)exp(-</2;tfcW).  (2) 


This  researen  was  supported  by  tne  Advroced  Researen  Protects 
Agency  of  the  Deoartmem  of  Defense  and  was  monnmad  by  me 
Air  Force  Office  of  Scientific  Researen  unoer  contract  No. 
F49620-90-C-0016.  The  United  States  Government  is 
autnonzeo  to  reoroouce  and  disznoote  reonnts  tor  governmental 
oumoaes  notwumtanding  any  cooyngm  notation  ncieon 


If  the  modulus  of  the  Zak  transform  of  w.  |Zwj,  is 
constant  almost  everywhere,  tfaei  wnnl  are  orthogonal 
and  well  behaved.  When  IZwj  has  a  small  dynamic  range 
(small  ratio  of  its  maximum  and  minimum  values),  the 
expansion  is  nonorthogonai  but  still  quite  stable.  The 
difficult  case  is  when  Zw  has  a  zero  within  0  £  i  <  T.  0  £  n 
<  UT.  The  formal  expression  for  the  Gabor  coefficients. 

XtTr.  Zf(  V.  T> 

a*j,  ■  I  f  —  expU2xinvT  -  mx  /  T)]dvdx,  (3) 

0  JZwfV.T) 

dons  not  than  necessarily  represent  the  Fourier  u  affluents 
of  an  function,  and  it  implies  an  expansion  for 

which  the  above  equation  does  not  always  yield  a  reliable 
basis  for  numeral  evaluation  of  theio^j. 

For  computation  by  the  Zak  method.  (3)  must  be 
replaced  by  a  finite,  discrete  equivalent.  Any  algorithm 
for  this  has  the  result 


Zfi  yr.) 

/m OgaO  vq.  “p) 


expC/2*tn  vfT  -  mxp  /  T)l  (4) 


Sampling  and  penodtzauon  can  be  applied  in  a  way  that 
mtaim  much  of  the  comispoooBitcc  to  the  cooumious  tunc 
expansion  fl.  2).  In  this  paper  we  consider  cransfonns 
made  discrete  by  sampling  and  truncation. 

Most  of  the  continuous  ante  accuracy/stability  issues 
carry  over  to  the  discrete  case.  In  particular,  it  is  easy  to 
see  in  (3)  the  impact  of  a  zero  of  Zw.  Even  if  Zf  and  Zw 
have  a  coinciding  zoo.  the  rauo  ZffLw  may  not  be  well 
defined.  When  tZwj  maxes  a  close  approach  to  zero  at 
some  grid  point  the  rauo  can  become  very  large  there, 
and  rtommain  the  expression  for  the  coefficients.  Ways  to 
control  transforms  with  Zak  zeros  are  given  in  secoon  2. 

A  second  approach  is  to  treat  the  sampled  equations 
for  the  Gabor  transform  as  a  matrix  and  solve  for  the 
coefficients  by  matrix  inversion.  The  structure  of  the 
matrix  is  such  that  the  inversion  can  be  performed 
efficiently  f3.  4).  For  reasons  of  completeness,  we  will 
now  introduce  the  rudiments  of  the  matrix  lormuiaoon.  If 
one  starts  with  the  Gaoor  transform  given  oy  CD. 
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uniformly  samples /«  every  &  i*T/M(tk»k(T/M)). 

and  truncates  the  expansion  after  P 
pouts,  one  ootains  the  represen  tao  on 

A1*)-  f  -*7)exp0‘2»«t  IT)  (5) 

aaOnO 

which  can  be  wnnen  in  matnx  form  as 

F-(W  E  )A  .  (6) 


If  Zw  has  a  zero  on  the  grid,  the  condition  moo  will  he 
infinite  and  the  transform  aoninveroblc.  Restricted  to  (9). 
C  need  not  be  infinite  when  (8)  is.  tet&ng  the  discrete 
transform  possibly  be  mote  stable  than  the 

Three  methods  of  stahiltnng  the  discrete  trensform 
am  illustrated.  We  soot  with  a  discrete  domain  winoow 
having  a  Zaic  transform  zero  on  the  gnd.  In  the  first 
method,  we  move  the  zero  of  the  Zak  transform  off  the 
grid  by  a  subtle  change  in  the  shape  of  the  window 
function.  The  second  method  is  to  apply  a  ume  and/or 
frequency  oansiaooo  to  the  anginal  wimow  that  relocates 
the  Zak  zero  according  to  the  formula 


where  F  and  A  are  P  x  1  column  vectors  corresponding  to 
the  signal  and  the  <  reordered)  coefficients  respectively, 
and  W  and  Eirefif  matrices  corresponding  to  the 
values  of  the  window  and  the  exponentials  associated  with 
the  rotations  around  the  unit  circle.  As  it  turns  out.  the 
matnx  W  is  block  lower  triangular  and  each  of  the  sub- 
blocks  is  diagonal.  The  sub-blocks  which  lie  on  die 
diagonal  are  all  equal  and  correspond  to  the  values  of  the 
window  over  the  first  M  points,  and  similarly  all  the  sub¬ 
blocks  on  the  f*  off-diagonal  are  equal  and  contain  the  M 
values  of  the  winoow  after  skipping  over  IM  -  1  points. 

The  matnx  £  also  has  a  particularly  nice  form.  It  is 
block  diagonal,  all  of  the  sub-blocks  are  equal,  and  the 
entnes  of  the  sub-blocks  are  the  entnes  of  the  Fourier 
rotation  matnx.  Consequently.  E  is  invertible  ani vocally. 
With  this  in  mud.  we  can  solve  (6)  for  the  unknown 
vector  of  Gabor  coefficients  and  we  obtain 

A-E^WF.  (7) 


z{/(r-t#)e>2~*‘'*t«>}- 
e>2«v.(r-r.>zJ|^  v_  v0,  r-r0).  (10) 

A  third  method  is  to  ream  the  window  function  but 
change  to  a  new  grid  on  which  the  zero  no  longer  lies. 

The  window  funcooo  employed  for  the  experiments  is 

w(r)  ■  rexp(— or  )u(r),  (11) 

where  a  is  a  positive  decay  constant  and  ins  the  umt  step 
function  at  the  angm.  The  Zak  transform  of  (11)  is 


Since  the  inverse  of  E  is  analytically  raimitwH  well  which  has  a  single  first-order  zero  at  the  pomt 
behaved,  the  stability  issues  are  solely  associated  with  the 

inverubiiiry  of  W.  The  structure  of  die  inverse  of  W  and  ,y-  r- .  _  f  1  T  \ 

the  resulting  stability  issues  will  be  discussed  in  section  3.  V  27' 7+7“^"  j 


(13) 


2.  ZAK  METHOD  EXPERIMENTS 
2J  Theory 

The  stability  of  a  continuous  tune  Gabor  representation 
can  be  gauged  by  the  srn  of  the  mno 

max  |Zw(  v.  r)!2 

^  OS  v<|/7\0Sr<r 

C  *  1  1  t  .  (8) 

noun  iZw<  v.  r)T 

osv<i/r.osr<r 

which  acts  as  a  condition  number  of  the  transform.  This 
number  is  also  the  frame  bound  ratio  of  the  transform  (5]. 

When  implemented  by  a  discrete  Zak  method,  the 
Gabor  trensform  is  characterized  by  the  same  ratio 
sampled  on  the  tone-frequency  gnd. 


Nonce  that  the  zero  lies  on  a  frequency  gnd  pomt  when 
the  number  of  frequency  points  M  is  even,  bat  the  tune 
location  of  the  zero  is  arbitrary  within  (772. 7). 

For  an  N  by  M  Gabor  trensform.  the  values  of  a  lying 
on  the  gnd  are  obtained  from  (13): 

u«yln^— —  lj;  0<m<Af/2.  (14) 

The  following  experiments  use  N  ■  16.  M  •  16.  Fig.  1 
shows  the  window,  and  Fig.  2  shows  the  Zak  trensform 
for  the  choice  a  «  -1.94591.  Because  of  smgie-precision 
roundoff  in  the  calculation,  the  condition  number  (8)  is 
8.523  x  1015  and  not  infinite  at  this  point  If  we  choose 
the  decay  constant  according  to  the  formula 


v«n/W7,r»«r /  Af:0Sn<N-l.0S/n<  M  —  \.  (9) 
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me  Zak  zero  will  lie  midway  between  grid  potnc.  For  tte 
cnoKX  a  ■  -1.68640.  we  find  the  Zak  tnumorm  shown  in 
Fig.  3  with  condition  numoerC  ■  4.175  x  I02 


Fig.  3  -  Zak  transform  of  the  analysts  wtnaow. 
era* re  the  analytical  zero  nas  oeen  moved  off  the 
grid  by  tul/urrmg  the  time  constant  The  condition 
member  u  displayed  on  trie  lower  left  hana  comer 


Fig.  i  ■  Analysis  wutaow 


Fig.  2  -  Zak  transform  of  the  window.  The 
transform  has  an  analytical  zero  on  a  time 
jrequencv  grid  pouit  The  condition  number  is 
displayed  on  the  lower  left  hana  comer. 

The  behavior  of  a  transform  having  a  Zak  zero  near 
the  grid  is  quite  sensitive  to  the  choice  of  a.  For  example, 
choosing  a  *  -1.94590.  a  change  in  the  5“  significant 
aigu.  changes  the  condition  number  by  4  oners  of 
magnitude  to  1.618  x  10"  When  the  zen  is  placed 
midway  between  grid  points  the  sensitivity  is  much  less, 
tor  example,  both  a  *  -1.68640  and  a  s  -1.68641  yield  C 
»  4. 175  x  103  to  four  significant  figures. 

Fig.  4  illustrates  one  appiicauon  of  our  second 
method:  the  window  having  the  on-grid  zen  is  tune- 
delayed  by  T/1M  tone-half  sample)  before  sampling,  and 
as  a  result  the  condition  number  (8)  of  the  Zak  transform 
becomes  C  » 1513  x  107. 

Experiments  showing  the  third  method  of  relocating 
the  zero  of  the  Zak  transtorm  ov  means  of  a  grid  change 
are  not  included  due  to  space  limitations 


Fig.  4  -  Zak  transform  of  tne  window  after  the 
analytical  zero  has  been  moved  off  the  gnd  bv 
delaying  the  window  arrival  tune.  The  condition 
number  u  snown  an  the  lower  left  hand  comer. 

3.  MATRIX  METHOD  EXPERIMENTS 
3.1  Theory 

As  we  mentioned  earlier,  the  stability  of  a  discrete  tune 
matrix  Gabor  representation  depends  only  on  tbe  inverse 
of  the  window  function  which  we  will  denote  by  B  due  to 
the  biorthogonaiity  relauonsnio  oetween  B  and  W  The 
form  of  B  is  identical  to  the  form  of  W.  and  its  entries  are 
given  by  the  recursive  relation 

B,  — W,-*(WfW0-'  -W,_IB1-r..^W1Bf_1)  (17) 

where  the  subscript  /  *  0. 1 . N-\  denotes  the  distance  of 

the  block  from  the  main  diagonal.  Gearty  the  only  part  of 
W  that  gets  inverted  is  the  block  associated  with  the  first 
M  points  of  the  winoow  function,  and  therefore,  a 
necessary  condition  for  mvernbiiity  is  that  the  window  be 
nonzero  everywnere  in  mat  region.  Even  if  that  condition 
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is  *  tt  does  not  imply  that  the  matrix  B  will  be 

well  beaaved.  Noace  that  as  we  proceed  funner  down  the 
off  diagonals.  ie.  as  L  gets  larger,  the  last  term  in  (17)  will 
be  multiplied  by  the  /*  power  of  the  inverse  of  W.  This 
essentially  means  that  comparatively  small  values  of  the 
window  over  the  first  M  points  can  lead  to  values  at  the 
tail  end  of  the  bioRhofooal  which  are  undesireably  large. 
A  measure  of  the  stability  of  the  window  can  at  least 
empirically  be  characterized  by  the  biorthogonal  having 
no  small  values  over  the  first  M  points,  compared  to  the 
values  over  the  test  of  the  window  values. 

12  Experiments 

The  window  function  used  for  the  matrix  experiments  is 
the  same  as  the  window  used  for  the  Zak  method 
experiments  with  the  time  constant  corresponding  to 
placing  the  zero  of  the  Zak  transform  on  one  of  the  grid 
points.  If  we  choose  this  window  as  is.  we  know  a-pnon 
that  the  biorthogonal  will  be  infinite,  since  our  window 
function  has  an  analytical  zero  at  the  origin,  and  this 
violates  our  initial  necessity  condition  for  invertibility. 
This  problem  can.  of  course,  be  alleviated.  If  we  shift  our 
window  to  the  left  by  2T/M.  ie  two  points,  the  window  no 
longer  lias  a  zero  on  the  interval  containing  the  first  M 
points,  and  the  biorthogonal  is  well  benaved  as  is  shown 
in  Fig.  5.  The  values  over  the  first  M  points  are  not 
significantly  smaller  than  the  rest  of  the  values,  and  the 
set  of  values  over  the  whole  interval  are  bounded  by,  in 
this  case.  13.0.  One  of  the  other  ways  to  obtain  a  nice 
biorthogonal  is  to  increase  die  amplitude  of  the  original 
window  by  a  constant  factor,  ie.  add  a  value  of  02  to 
every  point  of  the  window.  The  biorthogonal 
corresponding  to  this  window  has  values  that  are  bounded 
by  23.0.  and  it  is  displayed  in  Fig.  6. 


Fig.  5  •  Biorthogonal  to  the  analysis  window 
whose  zero  has  Been  shifted  two  tone  stunples  to  the 
left.  The  maximum  value  of  the  biorthogonal  is 
given  on  the  upper  left  hand  comer  of  the  was. 


alleviate  some  of  the  problems  associated  with  the  Zak 
transform  of  the  window  having  a  zero  somewhere  in  the 
region  of  interest.  As  we  have  shown  in  the  above 
experiments,  the  existence  of  a  second  method  to 
calculate  the  coefficients  does  not  guarantee  a  stable 
coefficient  set. 


Fig.  6  -  Biorthogonal  to  the  anatvsts  w tndow  eater 
adding  a  constant  value  of  0.2  to  everv  winaow 
sample.  The  maximum  value  of  the  btortnogotusl  Is 
given  on  the  upper  left  hand  comer  of  the  axis. 

One  has  to  resort  to  some  farther  knowledge  of  the 
structure  of  the  window  function  in  aider  to  oocun  results 
chat  are  not  sensitive  to  small  perturbations.  It  is  me 
belief  of  the  authors  that  by  appropriately  tailoring  the 
window,  the  instabilities  associated  with  one  or  both  of 
the  methods  can  be  minimized,  and  at  lean  one  of  the 
methods  will  give  satisfactory  results. 
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4.  SUMMARY 

The  initial  motivation  that  led  to  the  development  of  the 
matrix  representation  of  the  Gabor  transform  was  to 
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Abstract.1  Recent  interest  in  the  Gabor  transform  for  time-frequency  signal  analysis  can 
be  attributed  in  pan  to  increased  knowledge  about  the  accuracy,  stability,  and  complexity  of 
algorithms  for  computing  the  transforms.  Behavior  of  the  computations  depends  on.  among 
other  things,  the  manner  in  which  the  continuous-parameter  equations  are  made  discrete  and 
finite.  The  most  straightforward  means,  truncating  the  time  functions  to  compact  support  and 
sampling,  relinquishes  some  control  and  blurs  the  relationship  of  the  discrete  equations  to  the 
original  transforms.  A  more  satisfying  discretization  and  fimtizanon  process  that  preserves 
relations  to  the  continuous  parameter  case  is  found  by  periodization  and  sampling,  the  same 
method  used  to  obtain  the  finite,  discrete  Founer  transform  from  the  Fourier  integral.  By  this 
method  we  derive  the  finite,  discrete  Gabor  transform  equations  from  their  continuous 
parameter  counterparts,  in  the  process  explicitly  exhibiting  the  aliasings  that  permit  one 
periodic  sequence  to  be  the  finite,  discrete  Gabor  transform  of  the  other.  By  examining  the 
various  forms  in  which  the  Gabor  equations  can  be  expressed,  we  discover  how  the  input, 
window,  biorthogonai  function.  Gabor  coefficients  and  Zak  transforms  map  under 
periodization  and  sampling. 

1This  research  was  suooortad  by  me  Advanced  Research  Prefects  Agency  of  the  Deoartment  of 
Demise  and  was  monitored  by  the  Air  Force  Office  of  Scientific  Research  under  contract  No.  F49620- 
90-C-0016.  The  United  Stans  Government  is  authorized  to  reoroduce  and  distribute  rsormts  tor 
governmental  oumoses  notwithstanding  any  copyngnt  notation  hereon. 


1.  Introduction 


The  Gabor  transform  has  recently  attracted  a  certain  interest  in  time-frequency  signal 
processing  because  of  its  ability  to  portray  tune  history  of  the  frequency  content  of  a  signal  or 
image  by  representing  the  data  as  a  superposition  of  uniform  time  and  frequency  translates  of 
a  window  function.  Through  choice  of  the  window,  a  Gabor  representation  can  exploit  a  priori 
informauon  about  the  signal  that  is  useful  in  problems  of  signal  detection,  feature  extraction, 
and  identification. 

To  some  extent,  current  interest  in  Gabor  technique  is  fueled  by  successful  developments 
related  to  computing  the  transforms.  Since  Bastiaans  first  considered  the  question  of  soivmg 
for  the  Gabor  coefficients,  given  the  input  function  and  the  window  [9],  it  has  been  recognized 
that  this  computation  can  be  ill-posed,  leading  to  questions  of  accuracy  and  stability  of  the 
results.  Although  Bastiaans  illustrated  that  certain  choices  of  window  lead  to  particularly 
quirky  representations  and  provided  some  insight  into  the  mechanisms,  much  was  left  to  be 
categorized.  In  the  ten  years  following  his  valuable  work,  certain  of  the  practical  problems  in 
performing  the  computations  necessary  to  evaluate  the  transforms  have  come  under  careful 
scrutiny.  Introduction  of  the  Zak  transform  in  the  role  of  operational  calculus  for  the  Gabor 
theory  has  led  to  an  improved  understanding  of  the  underlying  mappings,  and  this  in  turn  has 
shed  light  upon  the  supporting  computational  processes.  Although  this  task  is  far  from 
complete,  much  more  is  now  known  about  the  accuracy  and  stability  of  algorithms  for 
computation  of  the  Gabor  transform  [5,  6].  In  addition,  the  computational  complexity  of  these 
algorithms  is  also  better  understood  [23,  24].  We  can  now  safely  pay  a  little  more  attention 
to  what  it  is  we  ought  :o  compute. 

In  deriving  digital  computation  algorithms  for  the  Gabor  transforms,  one  must  first  make 
the  equations  discrete  and  finite.  The  most  straightforward  means  is  to  truncate  the  time 
functions  to  some  compact  support  and.  sample  them  at  a  rate  that  captures  their  significant 
behavior.  This  process  relinquishes  some  control  over  the  end  product  and  blurs  the 


wiftny  pmmi. 

Nonetheless,  some  consistent  theones  can  be  developed  in  this  manner.  Wexkr  and  Raz 
[27]  have  suggested  a  discrete  formulation  in  which  the  ame  extent  of  the  window  function  is 
truncated  but  that  of  the  biorthogonal  function  is  not:  furthermore,  the  frequency  index 
remains  unbounded  Those  authors  have  commented  that  the  unbounded  supports  make  this 
form  ill-suited  to  digital  computation. 

More  recently,  Balart  [7]  has  developed  a  finite-dimension  matrix  formulation  of  the  fully 
truncated  case,  in  which  it  can  be  argued  that  the  "true"  discrete  Gabor  coefficients  are 
produced  when  the  window  has  bounded  support.  The  approach  of  Ausiander  et.  al.  [2.  61, 
which  uses  the  Zak  transform,  is  fully  discrete  and  finite,  and  produces  an  invertible  transform 
when  no  sample  of  the  Zak  transform  of  the  window  function  is  zero.  The  coefficients 
produced  by  this  algorithm,  however,  can  differ  considerably  from  those  produced  by  Balarf s 
method,  as  numerical  experiments  performed  by  Balart  and  the  author  have  shown.  Both 
methods  are  highly  useful  transforms  that  nevertheless  differ  slightly  from  the  continuous 
parameter  Gabor  transform.  Results  of  the  two  converge  for  windows  having  sharp 
discontinuities,  deviating  more  and  more  from  one  another  as  the  window  appears  to  become 
continuous.  This  is  not  unexpected  based  on  the  Zero  Theorem  of  Ausiander  and  Tolimien 
[4.  Theorem  IL2],  which  states  that  if  the  Zak  transform  of  a  sufficiently  rapidly  decaying 
function  is  continuous,  it  must  have  a  zero.  It  is  empirically  true  that  the  two  methods  are 
nearly  the  same  if  the  Zak  transform  of  the  window  has  a  narrow  dynamic  range,  but  they  can 
diverge  if  the  samples  of  that  Zak  transform  make  a  close  approach  to  a  zero. 

A  discretization  and  fimtization  process  in  which  some  of  the  relations  to  the  continuous 
parameter  case,  are  preserved  would  be  more  satisfying  and  would  lend  insight  into  how 
successful  the  more  brute-force  approach  might  prove.  Wexier  and  Raz  [27]  have  presented 
the  sampled  approach  for  periodic  functions— one  of  the  methods  discussed  in  this  paper — 
through  a  matrix  formulation,  and  Ausiander  et.  al.  [2,  6]  have  given  the  7-ak  transform 
relationships  under  periodization  of  the  signal.  The  contribution  of  the  present  paper  is  a 
unified  treatment  of  discretization  of  the  Gabor  transform  and  its  ancillary  funcuons  by 
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penodizing  and  sampling. 

The  Gabor  discredzaiion  issue  parallels  what  is  done  in  transitioning  from  the  Fourier 
integral  to  the  finite  DFT,  where  penodizing  and  sampling  in  the  time  domain  likewise 
and  samples  the  frequency  domain  picture.  Moreover,  the  two  periodic  sequences 
so  created  are  related  by  the  finite  discrete  Fourier  transform  (DFT).  The  result  for  a 
function  /  and  its  Fourier  transform  F 


F\v)  = 


dtfit)  exp(-)2jrw). 


is  that  the  following  two  sums  are  a  DFT  pair  [25]: 


(I  M  +  kP)T!P\\ 


j:  1F[ (q+jPm 


p,  qe  Z. 


(1) 


(2) 


For  this  rather  remarkable  formula — in  which  the  essence  of  the  Fourier  transform  property  is 
maintained  under  sampling  and  aliasing — to  bold,  the  sampling  rates  in  the  conjugate 
domains  must  be  properly  related.  One  samples  /  at  spacing  r,  =  T/P,  for  positive  T  and 
positive  integer  P,  and  the  resulting  sequence  is  aliased  at  every  F-th  sample.  The 
spectrum  is  then  sampled  at  interval  /,  *1/7  and  also  aliased  at  the  /*-th  sample,  resulting 
in  two  sequences  of  period  P.  The  relation 

Pfst,  =  1  (3) 

characterizes  the  parameter  constraint. 

It  is  important  to  keep  in  mind  that  equation  (2)  tells  us  nothing  about  the  resemblance  of 
these  finite  discrete  sequnces  to  their  parent,  continuous  parameter  functions;  however, 
conditions  under  which  one  is  essentially  a  sampled  version  of  the  other  are  well  known.  In  a 
clever  work.  Ausiander  and  Grtlnbaum  [3]  have  explored  signal-independent  error  bounds  for 
this  comparison,  in  the  process  numerically  verifying  the  necessity  of  condition  (3). 
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derive  ail  the  discrete  finite  Gabor  transform  equations  from  tbe  continuous  parameter 
versions  through  the  method  of  periodization  and  sampling,  in  the  process  explicitly  exhibiting 
the  aiiawngc  that  permit  one  periodic  sequence  to  be  the  finite,  discrete  Gabor  transform  of 
the  other.  By  camming  the  various  forms  in  which  the  Gabor  equations  can  be  expressed, 
we  discover  how  the  input,  window,  biorthogonal  function.  Gabor  coefficients  and  Zak 
transforms  map  under  periodization  and  sampling. 

Section  2  briefly  summarizes  tbe  continuous  parameter  Gabor  equations,  giving  the 
fundamental  definitions  for  and  relationships  among  the:  (i)  signal,  or  function  to  be  expanded, 
(ii)  window  function,  fiii)  biorthogonal  function,  fiv)  Gabor  coefficients,  (v)  Zak  transforms, 
(vi)  inner  products  in  the  time  and  Zak  domains,  and  (vii)  sampled  auto-  and  cross-ambiguity 
functions  of  the  signal  and  window.  In  Section  3,  we  transform  the  synthesis  equation  that 
recovers  a  function  from  its  Gabor  coefficients  into  discrete,  periodic  form.  The  reverse 
transformation  is  visited  in  Section  4  in  terms  of  the  biorthogonal  function.  In  Section  S,  the 
Gabor  relations  expressed  in  terms  of  Zak  transforms  are  discretized.  Behavior  of  the  inner 
product  of  Zak  transforms  under  the  discretization  process  is  shown  in  Section  6. 
Relationships  that  permit  solution  for  the  Gabor  coefficients  to  be  expressed  as  a 
deconvolution  of  the  sampled  crossambiguity  function  of  the  signal  and  window  are  analyzed 
■md-r  periodization  and  sampling  in  Section  7.  In  Section  8  we  obtain  the  Gabor  coefficients 
of  the  DFT  of  the  given  signal.  Our  concluding  observations  appear  in  Section  9. 
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2.  The  continuous  pars  motor  Gabor  eqnsd» 


In  this  section  we  summarize  without  proof  results  that  serve  as  a  basis  for  development 
of  the  discrete  equations. 

2*1.  Definition  (synthesis  formula) 

A  Gabor  representation  of  a  rime  function  f(t)  is  a  series  expansion  of  the  following  form: 


At)  *  X  Z  (4) 

mm -mm  n  mem 

in  which  fa*.*},  m,  n  e  Z,  are  the  Gabor  coefficients ,  and  w(t),  the  window. 
Translates  of  w(r)  over  the  von  Neumann  lattice  having  a  unit  area  cell  of  dimension  T  x 
(1  IT)  in  the  rime-frequency  plane  farm  the  Gabor  basis  { w^r)  ],  as  follows: 

wm.„(r)  a  w(t -  nTl  exp</2J0niD  \m*ne  Z.  (5) 

In  this  regard,  the  Gabor  expansion  may  be  compared  to  its  generalization,  the  Weyl- 
Heisenberg  expansions,  which  can  use  any  density  of  basis  functions  £  1  [13,  14,  18]. 

In  general,  the  {*,*„(*)}  are  nonorthogonal,  and  (4)  is  a  nonortbogonai  expansion, 
which  is  the  reason  one  must  take  some  pains  to  assure  accuracy  and  stability  of  algorithms 
that  compute  the  transformation.  Gabor  expansions  can  be  orthogonal,  but  this  is  neither 
required  nor  always  desirable,  as  discussed  in  [13].  Conditions  for  orthogonality  of  the 
Gabor  basis  are  found  in  Boon  and  Zak  [10].  These  can  be  expressed  in  terms  of  the  auto- 
ambiguity  function  of  the  window — see  Toiimien  and  Orr  [26] — and  a  general  construction  of 
orthogonal  Gabor  expansions  has  recently  been  developed  by  Coifxnan  ex.  al  [1 1]. 
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which  we  denote  by  g(t),  and  define  in  notation  different  from  Gabor's: 

M  _  (6) 

where  a  is  the  RMS  pulse  width.  The  normalization  in  (6)  is  that  the  window  has  unit' 
energy,  or  L\  norm,  which  is  a  convenient  condition  we  shall  assume  throughout  for  all 
windows: 

to  lg(*)P  *  1.  (7) 

Helstrom  [19],  in  1966,  demonstrated  the  continuous  parameter  version  of  Gabor's 
Gaussian  expansion  and  noted  the  relationship  to  Glauber's  coherent  states  of  quantum 
mechanics  [16].  In  1967  both  Montgomery  and  Reed  [22]  and  Crum  [12]  published  on  the 
expansion  of  functions  in  other  than  Gaussian  elementary  signals,  generating  what  was  then 
viewed  as  a  continuous  parameter  Gabor-like  transform  that  would  later  be  called  the  short- 
time  Fourier  transform. 

It  was  in  Basnaans  [8]  that  the  bi orthogonal  function  for  the  Gaussian  window  was 
introduced.  In  later  considering  other  functions  for  the  window  role,  Basnaans  [9]  provided 
the  basis  of  modem  Gabor  theory  by  his  concern  for  the  inversion  of  Gabor's  representation. 

2.2.  The  hiorthogonal  function  (analysis  formula ) 

In  the  hiorthogonal  method,  one  uses  both  the  basis  functions  { **.„( t) }  and  a  related  set 
of  hiorthogonal  functions  \  bm  „(t) ) . 

bmj ,(r)  =  bit  -  nT)  txpijltanjT) ;  m.  n  €  Z.  (8) 
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(wmji^p.ql  =  j  dt 


w(t  -  nT)bm(t  -  qT)  «p(/2JC(in  -  pWT]  a  6*^5^ . 


(9) 


When  such  b(t)  exists,  (9)  permits  analysis  of  /  into  its  coefficients  by  inner  products  of 
fit)  and  the  {£(*«(*))  according  to  the  formula 


Omji  *  {/>*■«)  *  dtj{t)  b^it). 


[ 


(10) 


2J.  The  Zak  transform 

A  second  method  of  writing  the  Gabor  coefficients  is  based  upon  using  the  Zak  transform,  a 
tune-frequency  mapping  given  by 


Zf(v,t)  a  £  f{kT  +•  t)  <xy(-j2izkVT). 


(ID 


Taking  Zak  transforms  on  both  sides  of  (4)  yields  the  following  relation  among  the  Zak 
transform  of  /,  the  Zak  transform  of  w,  and  the  Gabor  coefficients  (2,  6. 21]: 


Zf(v,  7)  a  Sv(v,t)  £  £  amjt  exp(f2ie(m^T- n\T)l 


(12) 


C-10 


The  {OmuJ  are  found  formally  by  mverunf  the  two-dimensional  Founer  sows  on  the  right 


Qm.n  3  J  drj  «p(/2*(-m*7V  nvT))-  (13) 

If.  however.  Zw  has  a  zero  where  Zf  does  not.  ZflZw  is  not  an  L2(TxF)  function,  where 
7  *  [0.  7)  and  F  *  [0,  1/7).  ZfTLw  can  be  Ll CTxF).  in  which  case  Gabor  coefficients 
can  be  computed  from  (13).  but  may  not  be  square  summabie.  and  any  truncated 
representation  that  uses  but  a  finite  number  of  them  can  have  large  errors.  Also,  a  Riemann 
sum  approximation  to  the  double  integral  may  not  provide  numerically  accurate  formulas. 
Issues  related  to  stabilization  of  this  process  are  not  discussed  in  this  paper. 

The  relationship  between  the  window  and  biorthogonal  functions  is  reciprocal  in  the  Zak 
domain:  Zb(v,  r)  =  l/T7Z’*w(v,  t )]. 

Synthesis  of  /  from  (11)  follows  by  inverting  the  DFT  after  recovering  Zflv,  t)  from  the 
[a^n)  Via  (12): 


r 

fikT  -r  r)  s  7  dv  Zfly,  x)  exp(/2jckvr),  k  e  Z.  re  T.  (14) 


2.4.  The  deconvoluaon  methods 

Deconvolution  permits  recovery  of  /  from  its  sampled  short-time  Fourier  transform  (STFT) 
{(/'  wm.n ))).  using  the  fact  that  the  Zak  transform  is  an  isomorphism,  preserving  inner 
products  to  within  a  scale  factor  [21]: 

WIZd=£*fl«). 


05) 


where 


£f! Zd*J  <*vj  dT  ^(v.OZ’gfv.T).  (16) 

In  [18]  and  [21],  the  authors  define  the  Zak  transform  with  a  normalization  that  makes  the 
scale  factor  in  (15)  unity  for  all  T.  Multiplying  both  sides  of  (12)  by  Z*w(v,  f). 


Z/(v.t)  Z’w(v.t)  *  |Zw(v,T)p  ^  X  a*jt  cxp(/2x(mt/T-  nvT)], 


(17) 


and  computing  the  Fourier  coefficients  of  both  expressions  reduces  (17),  alter  simplification, 
to  the  convolution  [2] 

(fl  =  X  X  0*> 

This  equation  can  of  course  be  obtained  more  directly  from  (4)  by  taking  the  inner  product 
indicated  on  the  left  and  employing  the  shift  theorem  [20].  For  present  purposes  it  is  the 
Fourier  transform  of  (17)  that  is  of  greater  interest. 

Equation  (18)  is  often  used  in  comparing  the  Gabor  and  short-time  Fourier  transforms 
[1].  It  shows  the  ‘smearing’  of  the  Gabor  coefficients  by  the  window-dependent  kernel  to 
produce  the  STFT.  If  a  function  w  is  used  as  a  window  in  both  transforms  it  cannot  have 
both  the  smoothness  desired  for  spectral  analysis  and  the  frame  structure  that  stabilizes  the 
Gabor  expansion,  as  explained  by  the  Baiian-Low  theorem  [17]. 

Extraction  of  the  Gabor  coefficients  from  (18)  requires  deconvolution  of  the  sampled 
STFT.  Fourier  analysis  of  both  sides  of  (18)  converts  (17)  into: 


C-  12 


Deconvolution  is  also  the  generalized  inverse  for  incomplete  Gabor  expansions,  providing 
the  coefficients  of  the  orthogonal  projection  of  /  onto  the  subspace  spanned  by  the  basis 
vectors. 

In  the  following  sections  we  convert  the  formulas  of  each  method  to  a  self-consistent  set 
of  equations  that  describe  the  finite,  discrete  Gabor  transform,  obtaining  tools  that  are 
amenable  to  digital  computation  of  the  transforms. 


3.  The  synthesis  formula 

In  this  and  the  following  sections  we  develop  the  finite  discrete  Gabor  relationships. 
Derivation  of  the  first  formula  is  presented  in  full.  Following  that  derivations  axe  reigaied  to 
the  Appendix. 

The  starting  point  is  (4),  from  which  we  eventually  want  to  prove  (33),  which  expresses 
the  samples  of  the  periodized  signal  (29)  as  a  function  of  suitably  penodized  window  samples 
(30)  and  Gabor  coefficients  (32).  We  periodize  f(t)  to  a  period  that  we  restrict  to  a 
multiple  of  T,  NT: 

m 

f/n(Os  'ZfU+kNT),  (20) 

*■  — 

and  formally  write  the  Gabor  expansion  of fmit)  as 

ffn(t)  =  XXX  4n./i  wit  -  nT  kNT)  exp  \J2Jon(t  -  IcNTYT).  (21) 

km-m  m  nm*m 
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If  we  make  the  following  rearrangement  of  (21), 


fm(Os 


1  I 


r  ■  n  i 


w[r-  (n  -  kN)T] 


iexp  \flTantfT], 


(22) 


we  see  that  cite  periodization  of  /  results  in  a  corresponding  periodization  of  window, 
wvr<f). 


wvriO*  X  w(r+  fcAT), 


(23) 


yielding 

«•  a* 

JyKO  =  X  Z  Bm*  "Mit-  nT)cxpij2nmt/T).  (24) 


By  (24),  the  Gabor  coefficients  are  invariant  to  this  joint  periodization  of  the  signal  and 
window. 

An  interpretation  of  more  immediate  use  results  from  a  rearrangement  of  (24).  Replace  n 
by  its  modulo  N  representation,  n  =  q  +  sN,  where  0  <,  q  <,N  -  l.jg  2,  and  replace  the 
corresponding  sum  by  a  double  sum  on  q  and  r: 

•  “  N-  1 

frd.t)  =  X  Z  X  "M<t-  qD^vUlnmt/T].  (25) 

mm~sw~m  <?»o 


The  periodicity  of  w  has  been  exploited  in  reaching  (25).  Reorder  the  summations  so  that 
the  innermost  is  over  s .  the  time  index  of  the  coefficients: 


/#r(0  *  X 


IAT10N  OP  THE  QABORIWPTOSEHT- 

^  N.  I  • 

*  X  X  am++*wta(t'  qT)t 

—  *■0 


AHON  TO  aflUTARY  PROBLEMS 


qT)  tx?\j2xmxfT[. 


The  bracketed  term  in  (26)  is  a  omc-aJiascd  Gabor  coefficient 


0$$  *  X  a*****  OZqZN-l. 


By  (26)  and  (27),  periodizingXf)  similarly  penodizes  the  window  and  the  {a*^}  in  time. 

Now  sample  ff/rit)  at  the  A/-th  harmonic  of  1/7*,  so  rh»r  the  jfc-th  sample  occurs  at 
time  kT/M.  The  resulting  sampled  sequence  then  has  period  P  -  MN,  and  (26)  becomes: 


jf  ■  X  i  «P  02 nmk/M) , 

m  m  -*•  <?■  o 


where  we  have  introduced  simpler  notation  for  the  pehodized  dam  and  window: 

AT  *  fmikT/M)  (29) 

=  WftrikT/Ml  (30) 

Now  apply  modular  representation  to  the  frequency  index  m,  since  the  complex 
exponential  in  (28)  has  period  Af  in  m.  Lei  m  =  p  *  rAf,  0£pSAf-l.  reZ,  and 
reorder  the  sums: 


ii"  =  I  I  S  exp  (P»Wp ♦  rMVim. 

/»  *  0  <7«  0  rM  — 
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Hus  operation  penodizes  the  Gabor  coefficients  in  frequency,  motivating  dm  definition 


X  ap™  rM.q  1  <1  I  *  riCq  -  sN  ■ 


(32) 


rnmhinmg  ail  the  above  we  have 


An  =  l'  t!  “fi,-*1  a&pxkptM) .  (33) 

pmO  <?m0 

In  summary,  the  finite  discrete  Gabor  representation  of  the  periodized  and  sampled  / 
uses  the  periodized  and  sampled  window  function  and  the  doubly  periodized  Gabor 
coefficients  of  the  parent  continuous  parameter  signal. 

4.  The  biorthogonai  formula 

This  method  begins  with  (10),  which  expresses  the  Gabor  coefficients  as  inner  products  off 
and  the  set  of  biorthogonai  functions,  and  executes  the  periodization  (32).  Following  steps 
similar  to  the  above  we  find 

olti”  *  (34) 

rm  0 

Thus  the  process  of  sampling  and  penodizing  replaces  the  integral  expression  for  a*.*  by  a 
finite,  discrete  inner  product  of  the  sampled,  periodized  data  and  a  similarly  sampled  and 
periodized  set  of  biorthogonai  functions.  The  factor  T/M  in  (34)  is  the  sampling  interval  and 
cakes  on  the  role  played  by  the  differential  dt  in  (10). 
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5.  tw»  y.PWLnBbfg^CA^nOW  OF  THE  GABOR  HgPHgaPfTATlOW  TO  MLfTARY  PROBLEMS 


In  (12),  which  states  the  relationship  of  the  Zak  transfonns  of  the  signal  and  window,  we  can 
replace  the  index  variables  as  follows: 

m  a  p  +  rM:  OZp&M  - 1,  rs,Z 

n  *  q  +  sN,  OZqZN- Use  Z  ■ 


and  sample  in  tune  and  frequency,  Le. 


v 


JL.  T  =  2lT 

NT'  M  ' 


(36) 


to  get 


Af  •  1  N  -  1 

OfU.rn  =  (Zwk.„X  X  »  ap^J2n(mp/M-nq/N)].  (37) 

P-0 q-0 

Above  we  have  introduced  a  shorthand  notation  for  the  Zak  transfonns. 


Zfln/NT.  mTIM)  m  (7/1^  .  (38) 

The  intermediate  steps  are  found  in  the  Appendix.  Inversion  of  the  2-0  Fourier  series 
immediately  yields 

<*’ *  *  i'  Hr*- -Vtfl*-**  -  nqm-  (39) 

p-Q  q- 0 

Next  we  show  how  to  recover  the  finite  discrete  version  of  a  function  from  its  Zak 
transform,  starting  from  (14)  and  periodizing.  The  periodization  effects  sampling  in  the 
frequency  variable  at  spacing  1  /NT.  Then  sampling  at  r » pT/M.  0  £  p  £  M  -  1,  yields 
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/«[(■£"►  *)rj  ft  ann,txp<l2titqm  0<p<M- Ik e  Z  (40) 

Formulas  (37)  and  (39)  show  that  the  discrete  Gabor  coefficients  will  exist  either  for  all 
/  in  LHZIP)  when  Zw  has  no  zero  on  the  sampling  grid,  or  in  general,  in  the  subspace  of 
functions  /  whose  Zak  transforms  vanish  at  the  zeros  of  Zw.  Via  (40).  the  P  samples  off 
can  be  recovered  from  the  P  samples  of  its  Zak  transform  by  a  set  of  DFTs.  Because  of  its 
quasi-periodicity,  the  Zak  transform  acquires  no  additional  periodicities  in  the  process. 

6.  The  inner  product  formulas 

To  streamline  the  derivation  of  the  finite  discrete  deconvolution  equations  in  Secuon  7,  we 
will  need  to  know  how  inner  products  map  under  periodization  and  sampling.  Let  and  g(P) 
be  two  penodized.  sampled  functions.  A  natural  definition  of  them  inner  product  is 


(/"  '*<">=  i  *  i'  tp  (*rr  (4i  > 

pmQ  p  *0 

We  find  that  the  result  may  be  written  in  either  of  the  forms 

(/r>ig<n)  =  ±Y  X  WWZfJL. 

.to 

-bi  1  («> 

m  *0  n m 0 

Thus  the  inner  product  of  two  penodized,  sampled  functions  is  proportional  to  the  2-D  inner 
product  of  their  Zak  transforms. 
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7.  The  SlRhMHMMkiCATION  OF  THU  GABOR  REPRESENTATION  TO  MKJTARY  PROBLEMS 


The  equation  for  the  deconvoluuon  of  the  pehodized.  ««wp ted  Gabor  coefficients  from  the 
sampled  short-time  Fourier  transform  is  readily  found  from  (33).  Take  the  indicated  inner 
product  and  apply  the  shift  theorem  to  the  indices,  yielding 

(/*!<?)*  X  X  aP%'N)  (wiri  1  w(*r-/>.s  -  }  •  (^3) 

pmO  qmO 

In  the  appendix  we  show  bow  this  result  may  be  found  via  the  Zak  transform. 

The  correspondence  of  (43)  to  (18)  is  clear.  The  doubly  infinite,  two-dimensional 
convolution  has  been  replaced  by  a  finite,  end- around  convolution.  Again,  the  doubly  aliased 
Gabor  coefficients  take  the  role  formerly  played  by  the  Gabor  coefficients  of  fit). 

8.  Gabor  transform  of  the  DFT 

The  Gabor  reprsentation  of  a  signal  and  its  Fourier  transform  are  tightly  interrelated  by 
way  of  the  Fourier  transform  of  the  window.  If  we  begin  with  the  discrete  representation  (33) 
and  take  the  DFT,  we  have 

fr*  *  X  s  X  I  X 

km  0  km  0  pm  0  tfmQ 

Performing  the  summation  on  k  first  isolates  the  expression 

N-  i 

=  £  wjP  cxpirj2icks/P),  (45) 

kmQ 

which  is  the  DFT  of  the  window.  In  terms  of  this  quantity,  the  transform  becomes 


<#{•*>  exp02it<mA(  -  r)k/P].  (44) 
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(46) 


fi*  •  X  X 

pm 0  faO 

In  (46),  the  r-th  DFT  coefficient  is  expressed  via  a  Gabor  expansion  whose  coefficients  are  - 
those  of/ and  using  window  w.  The  roles  of  time  and  frequency  have  been  interchanged  and 
there  is  a  sign  change  in  the  exponential.  This  formula  can  be  put  in  an  alternative  form  that 
recovers  some  of  the  lost  symmetry  by  replacing  q  with  N  -  n  and  denoting  the  Gabor 

a 

coefficients  of/  as 


On.  m 


JMJf) 
am.  AM' 


(47) 


Exrhangn  of  time  and  frequency  roles  is  now  complete;  the  coefficients  off  with  respect  to  w 
are  now  the  coefficients  of  /  rotated  90°  in  the  tune-frequency  plane: 


fr*  -  X  X  *  mNydN).  (48) 

naO *a0 


9.  Observations 

The  equations  defining  the  continuous  parameter  Gabor  transform  have  been  convened  to 
finite,  discrete  form  using  periodization  and  sampling.  The  periodic  sequences  {fi?)  and 
{<45*°).  where- 0  ZkZP  -  l,  OZpZM- l,  OZqZN  -  1,  become  a  Gabor  transform 
pair.  In  the  linear  transformation  relating  these,  the  window  function  is  replaced  by  the 
periodic  sequence  {w^},  0  £  k  £  P  -  1.  When  the  biorthogonai  function  is  well  behaved, 
the  f  aj,*’** }  can  be  obtained  from  an  inner  product  of  the  data,  and  a  sampled, 

pehodized  version  of  b.  Relationships  involving  the  Zak  transforms  map  similarly;  the 
continuous  arguments  are  replaced  by  discrete  ones,  Zf  — >  {(Z f)q,p),  without 
periodization.  Likewise,  the  convolution  equation  for  the  Gabor  coefficients  converts  to  a 


Fin*  THiAPPUCATlON  OF  THE  GABOR  REPRESENTATION  TO  MUTARY  PROBLEMS 
form  exhibiting  cross-ambiguity  functions  involving  the  periodic  {/*  }  and 

in  applying  the  Gabor  transform  to  signal  processing  problems,  one  would  like  to  have  a 

clear  correspondence  between  the  finite,  discrete  formulation  used  in  computation  and  the 

continuous  parameter  equations  of  the  theory.  Specifically,  one  wants  to  see  the  analog 

signal,  window  and  biorthogonal  functions  map  into  a  finite  set  of  their  samples,  and  to  have  a 

finite  set  of  Gabor  coefficients  that  (i)  are  defined  over  a  time-frequency  region  where  the 

coefficients  of  the  analog  signal  have  significant  values,  and  (ii)  are  approximately  equal  to 

the  latter  coefficients  in  that  region.  Creating  the  correspondence  by  periodization  and 

sampling  makes  this  relatively  straightforward.  Since  Gabor  coefficients  are  aliased  in  time 

and  frequency  in  the  discretization,  one  needs  to  assure  that  the  replicas  are  sufficiently 

spread  to  prevent  significant  overlap.  This  is  done  by  choosing  the  signal  period  to  be 

suitably  long  with  respect  to  the  lengths  of  the  signal  and  window  functions,  and  by  sampling 

at  a  rate  high  enough— -essentially  a  composite  Nyquist  rate — to  capture  the  significant 

frequency  behavior  of  the  signal  and  window.  A  sampling  theory  for  Gabor  expansions  will  be 

the  subject  of  a  forthcoming  paper  by  the  author,  but  in  its  absence  one  can  still  be  confident 

that  when  these  guideline  precautions  are  observed,  computations  made  with  the  finite. 

discrete  Gabor  transform  should  faithfully  reflea  the  behavior  of  their  continuous  parameter 

counterpart. 
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Z  bm  *  pkLnW  *  X  6w^Wexp(v2npiWDa  ^(f)  £  ex.p{-j2xpMtfT).  (A4) 

*>•-•  /»■  —  />■— 

The  summation  of  complex  exponentials  can  be  recognized  as  proportional  to  the  formal 
Fourier  series  of  a  train  of  impulses  separated  by  T/M.  The  r-th  Fourier  coefficient  of  such 
a  train  of  unit  impulses  is  given  by 
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Anal  TR  iAPPUCATION  OF  THE  GABOR  REPRESENTATION  TO  MUTARY  PROBLEMS 

cxpi-JlxrMtTT)  *  (A5) 


fVti 

a  « 

I  * 

X  Stt-kT7M) 

Jo 

**  — 

implying 

£  .  ^.w  -  *iu«  ($  i  *'-«»>■  <■ A6> 

/»--  *— 

When  this  is  into  (A3),  the  impulse  functions  sample  the  bioithogonal  funuons. 

yielding 

^  =  £1 (A7) 


where  is  the  sampled,  penodized  /,  see  definition  (29).  We  exploit  the  periodicity  by 
writing  the  index  k  in  terms  of  its  quotient  and  remainder  modulo  P ,  k  -  r  +  Ps, 
rearranging  the  sum  to  yield 


4S»  -  £  l‘  t 


(A8) 


The  bracketed  term  is  the  penodized.  sampled  biorthogonal  function  that  can  be  written  in  the 
notation  of  (29)  and  (30)  as 
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<A9) 


letting  us  write  the  final  result  as  in  (34). 


Beginning  with  (12),  and  making  the  change  of  variables  indicated  in  (35),  we  have 


M  - 1  N-  I 


Z/(v,T)  *  ZMy,x)  X  2  2  X  aP  +  *u,+  ,tfiXQ{j2x[fj>+  rMyt/T-  (q  +  sfif)vT)]}  (A10) 

rm-m  jm—  />■  0  7»C 


Sampling  (A10)  on  the  grid  given  in  (36)  yields 


M  -  1  N-  1 


X  X  X  X  (q  +  sAOn/AO] } *(A1 1) 

ra<o  ji«  *»0 


Upon  simplifying  the  exponential  terms  and  recognizing  the  inner  double  sum  of  the  Gabor 
coefficients  over  p  and  q  as  the  doubly  periodic  corfficients.  (A1 1)  becomes 

/ot<t+*D«  X  Ar*a-*NqT) 


*  ^  Tf 

JO 


Jv  9(v,  t)  exp(/2Jt(*  -  Aty)VT] 


*  7  dv  3f(v,  t)  exp(/2mt  vT)  X  V/H-  (A12) 


C-26 


Final  TRJfcPPUCAHON  OF  THE  GABOR  REPRESENTATION  TO  »BUT ARY"  PROBLEMS 

(%*(ZH *  (^k-^wL-  X  Z  "*W]-  (A21) 

F ■  0 


Now  sum  both  sides  of  (A21)  over  «  and  n.  The  left  and  right  sides  become,  respectively, 

Y  I  <A22> 

«  ■  0  u  ■  0 


and 


M-  1 


m  *0 


AT- I 

X  (2*4.*.  (Z*vr.*)IL*  X  X 

it  ■  0  />«0 


aj*^  OLp{j2x(mp/M-  nq/N)] 


s  X'  S'  aSS'M  X  X  (Zwk*  CZ»>r.  &j*  cxp(j2x(mp/M-  nq/N)] 

0m 0  <7*0  m  ■  0  n «0 


=  £  X  «iSJn  I  I  (a.u(zwr.^.jt- 

F*0  <7«0  ntaOnaO 


as  X  45”^  (W(fl  1  wMr.pj.) 

ft  mQ  qmQ 


(A23) 


Combining  the  above,  we  have  (43). 
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(A16) 


Attune  Mtwpiw  Bectrenica  Cocpoiatton 

N -  I 

2  exp{/2*(j  - «)  rW]  «  Vfi, . 

rmO 

so  that  (A15)  takes  the  form  given  in  (41). 

Using  the  easily  proved  relationship  [20] 

3f(v,fl  *  r/[-v,t)  (A17) 

and  the  periodicity  of  the  Zak  transform  in  its  frequency  variable.  The  Zak  transform  of  g  can 
be  expressed  in  terms  of  the  conjugate  of  g: 

(Zg)»J»»  3  (Zg  *  (Zg  W-A.MM  (A18) 


from  which  we  get  (42). 

Proofs  for  Section  T- 

We  muitipiy  both  sides  of  (37)  by  (ZwJL*,  finding 

V-IN-I 

(?A.m  (Z^  3  l(Z»^WP  Z  I  4^exp(/2w(mp/M-n^A0].  (A19) 

/7«  0 

If  we  muitipiy  (A19)  by  exp(/2jr(rp/Jlf  -  sq/N)],  the  left-hand  side  becomes 

WU  (ZH-t.  otpff2lt(rp/M  .  tqlN) ]  .  ( A20) 

where  the  second  factor  on  the  right  is  interpreted  as  the  (n,  m)-th  sample  of  the  Zak 
transform  of  wr<J .  The  iight-hand  side  of  (A19)  transforms  similarly,  leading  to 
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ABSTRACT 

Senes  expansions  of  signals  in  which  significant  features 
of  the  signal  are  captured  in  a  few  large  coefficients  are 
desirable.  This  work  shows  that  given  a  collection  of 
signals,  it  is  possible  to  find  Gabor  representations  for 
these  that  are  maximally  concentrated  in  tune-frequency 
space.  The  problem  addressed  is:  given  a  signal  set.  find 
the  window  function  of  the  Gabor  expansion  that 
minimizes  an  "average  dimension”  of  the  signal 
representations  relative  to  that  window.  The  dimension 
measure  employed  is  entropy  based  and  related  to  the 
quantum-mechanical  technique  where  one  interprets 
expansion  coefficients  as  probabilities. 

An  iterative  algorithm  based  on  partial  derivatives  of 
the  signal  set  dimension  with  respect  to  the  expansion 
function  was  used  to  evaluate  effectiveness  of  several 
nonlinear  optimization  algorithms  in  finding  an  optimum 
window. 


L  DIMENSIONALLY  OPTIMUM 
REPRESENTATION 

Senes  representations  of  signals  m  which  most  significant 
features  of  the  signal  are  captured  in  a  few  large 
coefficients  are  desirable  in  problems  of  detection,  feature 
extraction,  characterization  and  data  compression.  Many 
standard  expansions— «.g..  classical  Fourier  methods— 
have  but  little  flexibility  to  adapt  to  chanctensncs  of  the 
data.  Time-frequency  or  time-scale  representations 
encompassing  a  family  of  transforms  specified  by  a 
window  function  lor  analyzing  wavelet),  however,  have 
great  potential  to  accommodate  the  data  under  analysis  by 
selection  of  a  well-matched  transform.  To  date  there 
seems  to  have  been  little  effort  expended  in  achieving 
some  of  this  potential  for  economy  of  representation. 

A  theory  for  selection  of  good  basis  functions  is 
described  in  a  prior  publication  (1).  This  theory  centers 
on  defining  a  dimension— first  for  a  signal  and  then  for  a 
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signal  set— «n  terms  of  the  expansion  coefficients 
generated  by  the  signals  under  a  particular  representaoon. 
in  this  paper  we  report  early  experiments  in  implementing 
optimization  techniques  to  use  the  theory  in  obtaining 
ahnensi onaliy-opomum  Gabor  expansions. 

2.  THE  DIMENSION  OF  SIGNALS 
Time-bandwidth  product  is  a  familiar  measure  of  the 
complexity  of  a  signal.  The  BT  produo  is  intrinsic.  La.. 
independent  of  the  representation  applied.  By  contrast,  the 
measure  sought  here  is  to  be  used  in  preferentially 
selecting  among  representations,  and  thus  should  be 
representation  dependent.  The  same  should  be  true  of  the 
dimension  of  s  ret  of  signals,  although  such  a  measure 
begins  to  take  on  a  more  intrinsic  nature  as  suitable 
constraints  are  placed  upon  the  eligible  representations. 

Sines  dmcnanx  IC  hart  m  the  .»  frwrf 

all  its  properties  should  be  expressible  through  the  signal's 
expansion  coefficients.  When  the  expansion  is  a  Gabor 
representation  of  Sf(t), 
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where  the  (  w ^(t)  j  are  ame  and  frequency  translations 
of  a  window  function  w (r)  [2).  dimension  is  expressed  in 
terms  of  the  { a‘mjl  J. 

We  have  imposed  six  reasonableness'  requirements 
on  s  definition  of  dimension,  each  expressible  m  terms  of 
the  signal,  the  expansion  functions,  and  the  coefficients — 
see  [  1  ]  for  details.  Under  these  reauiremems  the  unique 
solution  for  the  dimension  has  the  form 
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The  result  is  an  entropy-derived  measure  for  which 
probabilities  are  developed  from  the  coefficients  I  Om.nh 
by  analogy  to  quantum  mechanical  methods.  When  a  set 
of  signals  {  Sf(t) )  is  considered,  a  set  dimension  can  be 
defined  as  a  weighted  average  of  individual  dimensions: 

t 

where  the  weights  {  di  I  could  represent  probabilities,  or 
acouid  be  assigned  to  reflect  relative  importance  of  the 
various  signals.  This  same  dimension  measure  has  been 
proposed  independently  by  Coifman  «  at  [3]. 


3.  SOLUTION  TECHNIQUE 
Analytic  minimization  of  (4)  is  in  general  a  bleak 
prospect.  Even  wriung  the  gradient  of  P  with  respect  to 
the  window  function  is  difficult,  and  the  solution  of 
VO  a  0  is  not  apparent.  We  address  the  minimization  as 
follows:  (1)  since  expansion  coefficients  are  linearly 
related  to  the  biorthogonal  function,  we  solve  not  for  the 
window  function  w  that  minimizes  D.  but  its  biorthogonal. 
b— permitting  the  anal  vac  gradient  to  be  obtained— from 
which  w  is  found:  (2)  solution  for  the  opamum  b  is 
accomplished  by  iterative  methods  of  nonlinear 
opumizauon  theory,  including  both  those  that  directly 
seek  the  minimum  of  D  and  those  that  solve  for  a  root  of 
the  gradient  of  P.  Derivative  aiding  can  range  from  none 
to  analytically  provided  second  derivatives  (Hessian 
matrix)  [4],  Routines  from  the  MATLAB  Optimization 
Toolbox  and  NAG  FORTRAN  Library  generated  the 
results. 

The  equations  resulting  horn  the  above  are 


d\nP 
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(5) 


To  do  the  implied  differentiation  of  the  (  pmjx }  requires  a 
relation  between  the  expansion  coefficients  and  the 
biorthogonal  function:  this  expression  differs  according  to 
the  form  of  discrete  Gabor  expansion  used.  Incorporation 
of  slope  information  in  seeking  the  zero  of  the  gradient  is 
achieved  using  the  mixed  second  order  paresis 


seven  iterations  to  the  almost  exactly  correct  state, 
achieving  a  dimension  of  1 .00001 .  Experiments  using 
multiple  signal  sets  and  other  algorithms  are  reported  in 
the  full  paper. 

Figures  2  and  3  illustrate  the  results  of  a  simple 
experiment  in  which  the  signal  set  contains  two  signals,  a 
rectangular  pulse  and  a  decaying  exponential.  Since  these 
two  have  distinct  envelopes,  no  one  wmoow  function  can 
represent  both  signals  with  unity  dimension.  Figure  2 
shows  the  representation  wnen  the  inonopumum) 
rectangular  signal  is  employed  as  the  window. 
Coefficients  for  the  exponential  signal  are  in  the  upper 
frame  and  those  for  the  puise  signal  in  the  lower  frame  in 
both  figures.  As  the  figure  snows,  the  Gabor 
representation  requires  3  large  and  several  small 
coefficients,  yielding  a  dimension  of  2.969.  The 
coefficients  for  the  pulse  signal,  as  expected,  consist  of 
only  a  single  term  ai  the  origin  with  a  dimension  of  1.0. 
The  average  dimension  for  this  representation  is  then 
1.98. 

The  window  found  by  opomizanon  (shown  in  Fig.  3) 
resembles  a  piecewise  linear  approximation  to  the 
exponential  signal,  and  the  biorthogonal  correspondingly 
resembles  that  of  an  exponential.  Gabor  coefficients  for 
the  exponential  signal  consist  of  one  large  term  at  the 
location  corresponding  to  the  starting  point  of  the  signal, 
and  several  small  coefficients  scattered  about  the 
frequency-tune  plane,  resulting  in  dimension  1.007.  The 
coefficient  map  for  the  pulse  signal  consists  of  two  large 
terms  at  the  first  two  Gabor  ome  pomts.  and  several  small 
coefficients,  producing  a  dimension  of  1.984.  The 
avenge  dimension  for  this  soiuaon  is  1J0.  which  is  25% 
lower  (better)  than  the  soiuaon  of  Fig.  L 

It  should  be  noted  that  we  have  not  shown  that  these 
results  represent  a  global  minimum  for  the  dimension.  It 
is  possible  that  a  more  sophisticated  cost  function  would 
produce  an  even  lower  dimension. 

Other  numerical  experiments  performed  with  various 
signal  sets  also  succeeded  in  finding  optima. 
Convergence  and  stability  of  soiuaons  were  found  to  be 
sensitive  functions  of  inihaiizauon.  signal  character! sues 
and  opdmizaaon  method.  Convergence  to  both  global 
and  local  minima  were  observed,  as  well  as  to  highly 
unstable,  asymptotically  global  minima. 


d2P  Ad\nPd\nP  ^  In D 
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4.  EXAMPLES 
Figure  I  shows  convergence  of  the  biorthogonal  function 
(solid  lines)  from  an  initial  to  a  final  state  when  the  signal 
set  consists  of  a  single  rectangular  pulse.  The  opamum 
-biorthogonal  function,  shown  dashed  in  the  initial'  panel, 
is  the  same  puise  in  this  case.  Starting  from  a  Gaussian 
random  vector  as  the  intuai  biorthogonal  funcaon.  the 
algorithm  (a  gradient-aided  steepest  descent  method  using 
mixed  polynomial  interpoiauon)  (4)  converges  within 
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